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ABSTRACT 

An extensive numerical experiment, using the developed computer code, has been 
conducted to design an optimized laser-sustained hydrogen plasma thruster. The plasma 
was sustained using. a 30 kW CO 2 laser beam operated at 10.6 fim focused inside the 
thruster. The adopted physical model considers two-dimensional compressible Navier- 
Stokes equations coupled with the laser power absorption process, geometric ray tracing 
for the laser beam, and the thermodynamically equilibrium (LTE) assumption for the 
plasma thermophysical and optical properties. 

A pressure based Navier-Stokes solver using body-fitted coordinate was used to cal- 
culate the laser-supported rocket flow which consists of both recirculating and transonic 
flow regions. The computer code was used to study the behavior of laser-sustained plasmas 
within a pipe over a wide range of forced convection and optical arrangements before it was 
applied to the thruster design, and these theoretical calculations agree well with existing 
experimental results. 

Several different throat size thrusters operated at 150 and 300 kPa chamber pressure 
were evaluated in the numerical experiment. It is found that the thruster performance 
(vacuum specific impulse) is highly dependent on the operating conditions, and that an 
adequately designed laser-supported thruster can have a specific impulse around 1500 
secs. The heat loading on the wall of the calculated thrusters were also estimated, and it 
is comparable to heat loading on the conventional chemical rocket. It was also found that 
the specific impulse of the calculated thrusters can be reduced by 200 secs due to the finite 
chemical reaction rate. 
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NOMENCLATURE 


A 

I 

J 

k 

P 

91 , 92,93 

9 

R, S,T 
s 

U,V 

u,t/ 

4 

p 

r 

p 

a 


coefficient in difference equation for variable 4> 
local laser ray intensity 

Jacobian of inverse coordinate transformation 

conductivity 

pressure 

geometric relations between coordinate systems 
optically thin radiation 

source term in equations 6, 7, 15, respectively 
distance along laser ray 

velocity component along f— and 77— axis, respectively 
velocity component along z— and r— axis, respectively 
axes of curvilinear coordinate system 
general dependent variable 
density 

effective diffusion coefficient 
viscosity 

absorption coefficient 


Subscript 
Ey W, N, S 
e,w,n,s 
i 

P 

rad 

<T,»7 

0 


four adjacent nodes to P 

four surfaces of the control volume centered at P 
for laser ray * 

nodal point to be solved in difference equation 
radiation 

partial derivative with respect to f , rj 
for dependent variable <f> 


Superscript 

* values based on given pressure field 

' corrected values according to corrected pressure field 
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I. INTRODUCTION 


A propulsion system supported by a laser-sustained plasma has been the subject of 
many studies. A schematic of a c.w. laser supported thruster is shown in Figure 1. The high 
power laser beam from a remote site is focused into the thruster. The plasma is initiated 
and continuous sustained near the laser focus. The plasma is used to absorb the power 
from the laser beam and convert it into the propellant. Since the focusing lens is positioned 
at the end of the thruster, the propellant flows into chamber radially through annular slots. 
Two slots are shown in the Figure 1. The propellant flowing through the slot near the lens 
is relatively cold and serves for lens cooling purposes. The major propellant enters the 
thruster through the other slot and has relatively high temperature. This temperature is 
dependent on the regenerative heating from the thruster wall to the propellant originally at 
the storage temperature. In the thruster, flow convects high enthalpy plasma downstream 
and mixes with other propellant to uniformize the radial temperature distribution. Then 
the propellant thermal energy is converted to the kinetic energy and produces the required 
thrust for the rocket through a converging and diverging nozzle. The major advantage of 
the laser thruster, compared to the chemical thruster, is higher specific impulse (about 
1000-2000 secs) since the lowest molecular weight propellant (hydrogen) and the higher 
propellant temperature (not limited to the flame temperature of a combustion process) 
can be used. 

The modeling of the laser-sustained plasma (LSP) for propulsion began with the one- 
dimensional work of Raizer. 1 Jeng and Keefer 2 had extensively reviewed the early analyti- 
cal work using one-dimensional and semi-two-dimensional models. A full two-dimensional 
model on the LSP was first numerically solved by Merkle and coworkers 3 ’ 4 . The objective 
of their work is to demonstrate the capability of the developed numerical method, and part 
of the governing process of LSP model were simplified. For example, the plasma radiation 
were neglected in their model. They succeeded in the calculation of a low temperature 
LSP (peak temperature about 4,000 K) in the hydrogen flow seeded with cesium. They 
adopted the rocket geometry; however, the calculations were limited to the subsonic flow 
region. It is feasible to extend their approach to a pure hydrogen plasma, but the higher 
peak temperature (about 16,000 K) and much stronger nonlinear optical properties of the 
pure hydrogen plasma can result in numerical stability problems. 

Jeng and Keefer 2 obtained the numerical solutions for a rigorous two-dimensional 
model with realistic gas properties for the LSP in pure argon subsonic pipe flow. Their 
results for the argon flows agree well with experimental data 6 ’ 6 in the prediction of plasma 
position, size, shape and peak temperature. Jeng and Keefer 7 ’ 8 also extend the model for 
the pure LSP in the pure hydrogen subsonic flow, and studied the effect of laser power, 
flow rate, static pressure and optical geometry on the LSP flows. A study using the 
different laser beam profile and wavelength of the potential developed high power laser 
were also investigated recently 9 . Their work indicated that there are no fatal errors in the 
laser propulsion concept, and the critical parameters, for example plasma size, position, 
radiation heat loss and energy conversion efficiency, can be controlled using an appropriate 
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'igure 1. Schematic of laser— sustained plasma thruster 








optical arrangement and a suitable flow configuration. 


The objectives of the current study are to extend the developed numerical code for 
calculating the realistic rocket thrusters, and to perform extensive numerical experiment 
for a optimized rocket design powered by a 30 kW CO 2 laser. In this report, the physical 
model and numerical algorithm used by the developed code and their verification will be 
first discussed followed by the thruster design study. 

n. THEORETICAL MODEL 


2.1 LASER BEAM 


The geometric optics was used to describe the laser beam which was assumed to be 
consisted from a finite number of individual laser rays. Each laser ray path through the 
optical system to the plasma flow was calculated using real-ray tracing. Diffraction of the 
laser beam due to the finite aperture of the lens, and effects due to the inhomogeneous 
refractive index within the plasma were neglected in the calculation. The local intensity 
of each individual ray t was described by Beer’s law written in the the following form. 


dU 

d&i 


= ali 


( 1 ) 


where s t - is the distance along the laser propagation direction for ray t, and a is absorption 
coefficient and depends on local plasma temperature and pressure. 

2.2 PLASMA FLOW 


The model assumed that the flow is laminar, steady-state, axisymmetric and has 
variable physical properties. The pressure of the plasma is relatively high so local thermo- 
dynamic equilibrium (LTE) for the plasma can be assumed with little error. Therefore, 
the plasma can be described by a single temperature, and its physical properties are only 
a function of this temperature and pressure. The model also considered the thermal radi- 
ation emitted from the plasma which included the contributions from free to free, free to 
bound, and bound to bound transitions. These radiation properties are a strong function 
of wavelength, and the current model adopted the simplified approach, developed by Kemp 
et al 10 and Caledonia et al. 11 to solve the hydrogen plasma radiation transport equations. 
Radiation transport was divided into two parts depending on wavelength: wavelength 
longer than 0.095 fi m, optical depth is long and an optically thin approximation, which 
all radiation escapes from the plasma, was used; and in other wavelength region, optical 
depth is short and an optically thick approximation, which can be described as a diffusion 
process, was assumed. Following the above assumptions, the equations of conservation of 
mass, momentum and energy for the flow can be written as : 
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where k e /f is the sum of molecular thermal conductivity and radiation-induced thermal 
conductivity in the optically thick limit; q ra d is radiative heat loss in the optically thin limit; 
aL Ii represents local power addition from laser beam to the flow; and as mentioned before, 
the physical properties involved in equations (2-5) are a function of the local temperatures 
and pressures, and Ref. 7 has a detailed description of the thermophysical and optical 
properties sources used in this work. 
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In the thruster calculations, constant wall temperature was assumed. Total pressure, 
total temperature and flow streamline were specified at the inlet, and upstream running 
Riemann invariants extrapolation plus the isentropic relations were used to calculate other 
properties. At the exit, the properties were linearly extrapolated from the interior flow 
field. 

m. NUMERICAL ALGORITHM 

An extensive discussion of the adopted numerical algorithm can be found in Appendix 
(A). Only a brief description of this algorithm is presented in this section. 

3.1 COORDINATE TRANSFORMATION 

The set of transport equations (3-5) can be written in the cylindrical system coordi- 
nates for the dependent variable, <f>, in the following generalized form: 

±(rpn<f,) + A(r„*) = 

+R^(x,r) (6) 

where T* t ff is the effective diffusion coefficient and R^(x,r) is source term. The above 
source term for the momentum equations (3,4) also includes the diffusion terms from the 
dilation effect. When new independent variables f and 17 are introduced, Eq. (6) changes 
according to the general transformation f and 77. A schematic illustration of the relations 
between the physical domain and the transformed domain is shown in Figure 2. Equation 
(6) is rewritten in f and 17 coordinates as follows: 

1 d . 1 d . Trj . 
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1 d rr +r, . . ,1 Id \T*r, 

Jdl [—(«»*« _ + q3 *"\ 


+s*((,<l) (7) 

where 

U = ur„ - Vi,, (8) 

V = vx^ — ur( (9) 

9i = A + rj (10) 

92 = XiX n + (11) 
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( 12 ) 

( 13 ) 


93 = x\ + y| 

J = x € r„ - x„r£ 

and <S"^(£, 77) is the source term in the f and 77 coordinates. 
3.2 FINITE DIFFERENCE METHOD 


A nonstaggered grid system, as shown in Figure 2, is adopted in the calculations. For 
a typical grid node P, where u, v,p and h five unknown are located, are enclosed in its cell 
and surrounded by its neighbors N, S, E, and W. The integral of Eq. 7 over its control 
volume, by arbitrarily taking Af = 1 and Arj = 1, becomes 


{rpU<j>) e - {rpU<j>) w + ( 1 rpV<j>) n - {rpV<f>) 8 


rT+ 

= “j“(9i0£ ~ 92^) 


rT+ 
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>r * l ivr* 




( 14 ) 


The second-order center differencing is used to discretize the above equation except for 
the convective terms when cell Reynolds number is greater than 2. In that case, the first 
order upwind differencing is used. A staggered grid system in the body-fitted coordinates 
has been widely used in the literature, 12 and it has been known that the oscillatory field 
will be produced using the nonstaggered grid system. The artificial viscosity term proposed 
by Rhie and Chow 13 to cure the problem was also implemented in this study. The cross 
derivatives in the diffusion terms of Eq. ( 14 ) are the results of nonorthogonal coordinate 
system and are usually small. In order to use the solution procedure for the five-point 
unknowns, these terms are combined with the source term and treated as known quantities 
and the resulting relation between <j>p and the neighboring values can be written as : 


A p <j>p = Ae4>e + An<!>n + As<i>s + S*r • T*(f, rj) J, (15) 

where the coefficients A involve the flow properties of convection, diffusion, area, etc, and 
T^(f, 77) includes the source term S^(f, 77) and the cross derivatives in the diffusion terms. 

3.3 PRESSURE CORRECTION PROCEDURE 


While the above general transport equation, Eq. ( 15 ), was used to solved momentum 
and energy equations, the pressure field must be solved using the continuity equation. The 
procedure used in this work is an adoption of algorithm developed by Rhie, 14 and only a 
brief discussion is presented in this section. The pressure equation used in the adopted 
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procedure is similar to that of SIMPLER 15 and PISO, 16,17 except the density term is 
treated implicitly for the compressible flow. The basic idea is to use the following relations 
to correct the mass flux obtained from momentum equations. 

pu = p*u* 4- p'u * + p*u ' (16) 

pv = pV+/>V-*-pV (17) 

and the discretized pressure equation becomes a transport equation for compressible flow 
instead of Poission equation. The ratio of diffusion term to convection term is dependent 
on Mach number, and as Mach number increases, the convection term becomes dominant. 
The split operator concept similar to PISO were also adopted in the current study to solve 
the pressure correction equation. 

IV. MODEL VERIFICATION 

Because of the lack of hydrogen plasma data, the model evaluation procedure is based 
on the existing argon plasma data. The model has been compared with a wide variety of 
experimental results given by Welle et al 5 - 6 in order to learn the capability and limitations 
of the model. Both calculations and experiments were limited to the LSP within a constant 
area pipe. Detailed discussion on the calculation procedure can be found in Refs. 2 and 7. 
Since the spectroscopic measurements were limited to temperatures greater than 10,000 
K, the following figures contain only those isotherms within that temperature range, and a 
limited portion of the calculation domain is plotted in order to gain a clearer comparison. 
Only a few typical comparisons are presented in this section, and extensive discussion of 
the model performance can be found in Refs. 2 and 7. 

4.1 INLET VELOCITY EFFECT 


The performance of the model for two different inlet velocities, using an 8-in. focal 
length lens, is demonstrated in Figures 3 and 4, respectively. The static pressure and laser 
power of these two LSP’s differ by less than 1 percent, and the inlet velocity of the LSP in 
Figure 3 is about seven times that in Figure 4. The model calculates the size and position 
of both plasmas well. As the inlet velocity increases, the radius of the plasma decreases, 
and the plasma is located further downstream. Two local maximum temperatures within 
the low-velocity LSP are predicted by the model. One is on the laser path upstream of 
the focal zone, and another is at the focal zone. Both calculated maximum temperatures 
are about 1,000 K higher than the measured peak temperature (14,500 K). Only a single 
maximum temperature was measured, and it is located neither on the laser path nor at the 
focal zone. The predicted local peak temperature at the focal zone, which was not observed 
in the experiment, may result from the neglect of diffraction and refraction of the laser 
beam. Although it has been demonstrated experimentally that this assumption leads to 
only a small effect on the measured global laser power absorption, the predictive model can 
suffer more serious effects. If refraction of the laser beam within the LSP is considered, it 
results in a larger focal spot size and reduced laser intensity. The spot size for the refracted 


8 




o 



(lulu) 83UD|S|p | d i p dj 


i n 


Figure 4. Temperature contour plot of a low-velocity plasma (contour lines 
starting from 10,500 K with 500 K increments) 



beam is highly dependent on the LSP position and size. For a low-velocity LSP, most of 
the plasma is located upstream of the focal zone, and this makes the focal spot larger. 
When diffraction is also considered, the diffraction limited focal spot is approximately 
twice the size of that calculated from a geometric ray trace. The combination of laser 
beam refraction and diffraction, and the longer absorption pathlength ahead of the focal 
zone results in an actual laser intensity which is smaller near the laser focal zone than 
that which is predicted by the model. Thus, the local temperature maximum will not be 
generated there in the low-velocity LSP. This may explain why the calculations are not 
consistent with the experiments in this detail. 

Since the size and temperature of these LSP were well predicted, it is not surprising 
that the calculated thermal radiation power loss (316 and 445 W for high and low veloc- 
ity, respectively) agreed with the experimental results (326 and 392 W for high and low 
velocity, respectively). The model prediction for the transmitted laser power in the high- 
velocity LSP also agrees reasonably well (110 W and 159 W for prediction and experiment, 
respectively). However, the predicted transmitted laser power (91 W) in the low- velocity 
LSP is much smaller than the measurement (246 W). This is because the predicted second 
local maximum temperature near the focal zone absorbs an unrealistically large amount 
of additional laser power. 

4.2 LENS EFFECT 


Isotherms are shown in Figure 5 for an LSP with operating conditions similar to 
those of Figure 4 except for use of a 12-in. focal length lens instead of the 8-in. lens. 
Both calculation and experiment show that the plasma moves further upstream than in 
Figure 4, and the model predicted the upstream position and the radius of the plasma 
well. Again, the model predicts a locatl maximum temperature near the focal zone, which 
was not measured in the experiment. The effect of laser beam refraction and diffraction, 
which was discussed in the last section, is probably responsible. Similar to the 8-in. LSP, 
the thermal radiation from the LSP is well predicted, and the transmitted laser power is 
underestimated. 

4.3 LASER POWER EFFECT 


Figures 6 and 7 show the effect of laser power on the structure of the plasma. Both 
cases use the 12-in. focal length lens with similau: flow conditions. The calculations aire 
quite consistent with the experiments. As laser power increases, the plasma becomes lairger 
and moves upstream. For the 10,500 K isotherms, the model underpredicted the length. 
However, the experimental results exhibit more noise in these isotherms than in the higher 
temperature isotherms, atnd if the 11,500 K isotherms axe used then the model predicts the 
length of this isotherm within 10 percent. The calculated positions of both plasmas aire 
about 1 mm downstream of the experiments. After carefully reviewing the experimental 
data, it was found tb»t the on-axis maximum temperatures are not on the laser beam 
path (using real ray tracing and neglecting the inhomogeneous refractive index effect). If 
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■ Figure 7 • Temperature contour plot of a plasma with low laser power (contour 
lines starting from 10,500 K with 500 K increments) 




the experimental results are shifted downstream by 1 mm, which is within the accuracy 
claimed by Welle et al, the locations of peak temperatures are on the laser path. In this 
case, the model predictions agree well with the experimental LSP positions. The predicted 
global properties are also in reasonable agreement with the experiments: the transmitted 
laser power is within 10 percent, and the thermal radiation from the LSP is underpredicted 
within 30 percent. 

4.4 STATIC PRESSURE EFFECT 


All the LSP’s shown in the previous sections have either 2.0 or 2.5 atm static pressure, 
and the model performed well at these pressures. In this section, the effects of the static 
pressure are presented. The intrinsic properties of plasmas are highly pressure dependent, 
especially the optical properties. The absorption coefficient of the plasma at the laser 
wavelength and the thermal radiation from the plasma are approximately proportional to 
p 1-5 . Figures 8 and 9 show the plasmas at 3.0 and 1.5 atm, respectively. Both plasmas 
have similar laser power levels and inlet mass flow rates (but different inlet velocity). In 
the low-pressure case, the predicted plasma size is much smaller than the experiment. In 
the 3 atm case, the model provides an excellent prediction. 

After reviewing the model assumptions, we concluded that the treatment of the opti- 
cally thick thermal radiation was the weakest link, and was probably responsible for the 
failure of the model in the low-pressure case. As mentioned before, the validity of the 
diffusion approximation for short range radiation is dependent on pressure and the size of 
the plasma. As the pressure decreases, the required size of the plasma increases. For the 
present calculation, the size of the 1.5 atm plasma is below the minimum size required for 
valid use of the diffusion approximation found by Kopainskyl8. Although the diffusion 
approximation for short-range thermal radiation is not valid for the low pressure plasma, 
a sensitivity study using this approximation was performed to investigate the sensitivity of 
the solutions to radiative transport. For small, low-pressure plasmas, the radiation induced 
thermal conductivity should be smaller than that for a large, low-pressure plasma. 

A calculated result, where the radiation induced thermal conductivity was reduced 
to 50 percent of the value used for the calculation in Figure 8, is illustrated in Figure 10. 
The size of this plasma is much larger than the previous calculation, and agrees well with 
the experimental results. The role of a change in thermal conductivity on the results can 
be understood qualitatively. As conductivity decreases, the conduction heat transfer loss 
from the portion of the plasma along the laser path is reduced, which increases the plasma 
temperature along the laser path. Due to this increased temperature, the absorption 
coefficient at the laser wavelength increases and the plasma absorbs more laser power, 
which in turn generates a larger plasma. The exact solution of radiation heat transport will 
not be easy to achieve in the near future, and the optically thin and thick approximations 
for the radiative transport process should be retained for the near term. Although in 
some pressure and dimension ranges, the optically thick approximation is not quite valid, 
this concept can still be utilized with some adjustment for the magnitude of the radiation 
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Figure 9 . Temperature contour plot of a. high-pressure ( 3.0 atm) plasma (contour 
lines starting from 10,500 K with 500 K increments) 
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Figure 10. Temperature contour plot of a low-pressure (1.5 atm) plasma using 

50 $ of radiation induced thermal conductivity (contour lines starting 
from 10,500 K with 500 K increments) 




induced thermal conductivity. 


V. LASER-SUSTAINED HYDROGEN PLASMA 

The behavior of laser-sustained hydrogen plasmas within a constant area pipe were 
also extensively studied. Parametric study includes the use of different wavelength lasers 
as well as various forced convections. Due to the lack of experimental data, only theoretical 
calculations are presented. The detailed structure of LSP using either hydrogen or argon 
are similar, and detailed discussion on this aspect can be found in Ref. 19. In the following, 
only the global energy conversion efficiency from laser beam to the flows and the LSP 
structure using different wavelength lasers are discussed. 

5.1 POWER CONVERSION EFFICIENCY 


Figures 11, 12 and 13 show power conversion efficiencies under various forced convec- 
tion conditions. At 1 atm a large amount of laser power is transmitted through the LSP. 
Transmitted laser power, which decreases as pressure increases, is less than 2 percent in 
the plasmas with a static pressure greater than 3 atm. The thermal power radiated in 
the constant inlet velocity cases decreases as pressure increases, but the thermal power 
radiated first increases with pressure, then decreases for pressures greater than 5 atm in 
the constant mass flow rate cases. Generally speaking, the energy conversion efficiency 
has only a weak dependence on pressure in the constant mass flow rate cases. However, 
in the constant velocity cases, the energy conversion efficiency is greater at higher static 
pressure. Influence of laser power and force convection on laser power conversion efficiency 
is shown in Figs. 12 and 13, respectively. If the laser power is larger than 10 kW and 
the static pressure is higher than 3 atm, the LSP absorbs nearly all of laser power. In 
these situations, the optically thin radiation heat loss is the critical factor in the energy 
conversion efficiency. 

5.2 WAVELENGTH DEPENDENCE 


All of the known experiments with continuous laser sustained plasmas have used the 
10.6 micrometer wavelength from the carbon dioxide laser. However, it now appears that 
lasers capable of providing the megawatt powers required for a practical propulsion system 
are likely to operate at shorter wavelengths. The high power lasers which are currently 
under development and are candidates for laser propulsion are either free-electron lasers 
or chemical lasers. Free-electron lasers can be operated over a wide range of wavelengths 
and could be tuned to the atmospherical transmission window near 2.2 micrometers for 
a ground based system. The chemical lasers operate in the wavelength range from 3.3 to 

4.3 micrometers, depending on the particular chemical reacting system. We have used our 
computer model to predict the characteristics of plasmas sustained with three different 
wavelengths: including 10.6 micrometers as a benchmark case, and 3.9 micrometers and 
2.2 micrometers as representative of the free-electron laser and chemical laser, respectively. 
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Figure 11. Fractions of the incident laser power that are absorbed 
and radiated for several static pressures 
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Figure 12. Fractions of the incident laser power that are absorbed and 
radiated for a range of incident laser powers • 
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Figure 13 • Fractions of the incident laser power that are absorbed 
and radiated for different inlet velocities 
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Figure 14. Temperature contour plot of plasma using 
three wavelength lasers 
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The major interaction between the laser beam and the plasma is the absorption of 
laser radiation by the plasma. For the three laser wavelengths considered in this study, 
plasma absorption occurs through an inverse bremsstrahlung process consisting of electron- 
ion, electron-atom, and ion-ion absorption. The absorption coefficient of the plasma in the 
LSP temperature range (less than 20,000 K) is approximately inversely proportional to 
the square of the wavelength. This relationship and the absorption coefficient at 10.6 
micrometers adopted in Reference 8 were used in the cuirent study. 


Calculated isotherms for three plasmas are shown in Figure 14(a-c). The incident laser 
power for the different wavelengths was scaled inversely proportional to the absorption 
coefficient. The incident laser beam profile and the radius of the lens were the same as 
those used in the Ref. 8, but the incoming laser beam diameter was chosen to be five times 
larger to keep the power density at the lens within reasonable values. This results in an 
f/number of 3.46 which was constant for each of the three cases. 

The scaling factor for the incoming laser power that was used for this study (in- 
versely proportional to the square of the wavelength), was chosen so the local laser power 
absorbed by the plasma (product of local laser intensity and absorption coefficient) was 
approximately the same for the three wavelengths. The attenuation of the laser beam 
intensity depends on the product of the absorption coefficient and the physical length. For 
the same propagation distance through the plasma, a greater fraction of the incident laser 
power will be absorbed for the longer wavelength laser. This difference in the character- 
istic absorption length is reflected in the inverse relationship between plasma length and 
wavelength shown in Figure 12(a-c). The extended tail of the plasma downstream of the 
focus is sustained by the larger fraction of incoming laser power that passes beyond the 
focus for the shorter wavelengths. The fraction of the incident laser power that is absorbed 
by the plasmas was highest for 10.6 micrometers among the three test cases, but the total 
laser power absorbed by the plasma was the smallest. 

VI. ROCKET DESIGN AND PERFORMANCE 

The extensive numerical calculations have been performed to have optimized thrusters. 
The preliminary rocket design effort using the developed code are reported in Ref. 19. 
In the following section, the discussion will limited to four thruster designs. The test 
conditions and geometry of theses thrusters are shown in Figure 15 and Table 1. The 
calculated isotherms and normalized streamlines for rocket No. 1 are shown in Figure 
16 and Figure 17, respectively. The plasma is positioned near the throat with a peak 
temperature of 15,800 K upstream of the laser focus. Most of the propellant within the 
thruster is not directly heated by the laser beam. However, nearly all of the propellant 
mixes with the hot plasma and raises the propellant bulk temperature up to approximately 
10,000 K before it enters the diverging section of the thruster. Eighty four percent of laser 
power (25 kW) is absorbed by the propellant, and most of the remaining laser power is 
transmitted through the throat, except for a very small portion that reaches the thruster 
wall. 
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Figure 16. Isotherms for rocket No. 1 (starting from 2,000 K to 15,000 K 
with 1,000 K increments) 




Figure 17. Normalized streamlines for rocket No. 1 (starting from - 0.09 
to 1,01 with 0.1 increments) 


Since the propellant enters the chamber radially through a slot away from the lens, 
part of the propellant is expected to recirculate within the thruster. As shown in Figure 
17, two recirculation zones occur in this thruster. A large one exists near the focusing lens, 
and the other one is near the thruster wall downstream of the inlet slot. The streamlines 
near the centerline upstream of the plasma indicate the flow has an outward direction. 
This is because the downstream hot plasma propellant (low density) generates a locally 
unfavorable axial pressure gradient in the flow. This phenomena was also observed in the 
previous work 2 ’ 7 ’ 8 for the LSP within a constant area pipe. 

An enlarged view of the isotherms and constant axial velocity contours around the 
throat are shown in Figure 18 and Figure 19, respectively. The leading edge of the LSP is 
11 mm upstream of the throat. The axial temperature gradient at the leading edge is very 
high since the axial heat transfer of convection and conduction is in opposite directions. 
The detailed discussion of how LSP structures are influenced by forced convection and 
optical geometry can be found in the previous work. 2 ’ 7 ’ 8 At the throat the propellant is 
very hot, approximately 13,500 K on the centerline. The temperature decreases radially 
from the centerline and reachs a value near 1250 K, which corresponds to the sonic tem- 
perature of the adiabatic (without heat addition from the laser beam) propellant. Then 
the temperature increases due to the wall heat transfer and viscous dissipation effects. 

The traditional chemical rockets are preferred to have a uniform temperature distri- 
bution at the throat to minimize the two-dimensional loss of thruster performance (Isp). 
For rockets supported by a laser beam, the spatial distribution of heat addition is highly 
nonuniform, and a nonuniform temperature distribution is likely to exit at the throat as 
shown in Figure 18. Actually, it is required to have some kind of nonuniform temperature 
distribution around the throat for a practical high performance laser supported thruster 
(around 10,000 K propellant temperature at the throat), or the thruster wall will be dam- 
aged. The heat loading on the thruster wall, which will be discussed later in this paper, is 
primary due to irradiation from the plasma. A stream of cold propellant near the wall as 
shown in Figure 18, which results in some loss of rocket performance, is required to cool 
the wall. 

Upstream of the leading edge of the plasma a local minimu m in the axial velocity 
occurs on the centerline where the propellant is unheated by the laser beam. The radial 
distribution of the axial velocity within the center plasma core region is relatively uniform 
with the peak axial velocity on the centerline. Large radial variations of the axial velocity 
is found near the thruster wall where the large radial temperature gradient is also located. 
The exit temperature and axial velocity on the centerline are 10,700 K and 25,300 m/s, 
respectively. The predicted vacuum specific impulse for this rocket, with exit area expan- 
sion ratio of 4, is 960 secs. With further expansion of this rocket, the estimated specific 
impulse of this rocket can be approximately 1500 secs. 

The isothei ms for two other throat radius thrusters operated at 150 kPa pressure, 
rockets Nos. 2 and 3, are shown in Figure 20. The leading edge and size of these plasmas 
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Figure 18. Enlarged view of isotherms for rocket No. 1 (starting from 2,000 K 
to 15,000 K with 1,000 K increments) 
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Figure 20-1. Isotherms for rocket (a) rocket No. 2 (b) rocket No. 3 (starting 
from 2,000 K to 15,000 K with 1,000 K increments) 
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2,000 K to 15,000 K with 1,000 K increments) 



are different. The mass flowrate, which is directly related to the throat size, is the major 
controlling factor for the plasma size and position. The mass flow rate of a constant total 
enthalpy fluid through a choked nozzle is proportional to the throat area. For the laser 
supported thruster, the propellant temperature and its distribution at the throat are other 
controlling factors for the mass flow rate. For a smaller throat thruster, a larger portion of 
the propellant at the throat has a high temperature and this will further reduce its critical 
mass flowrate. The calculated mass flow rate for cases 1-3, tabulated in Table 2, result from 
this effect. For a smaller throat thruster, the reduced axial velocity within the subsonic 
converging section of the thruster allows the plasma to move further upstream along the 
laser beam and to obtain a larger size, as shown in Figure 20. A detailed discussion of 
the interaction between the plasma structure and a forced convection flow can be found in 
References 2,7 and 8. 


In the traditional chemical rocket, the total pressure has a minor effect on the ex- 
pected value of specific impulse. However, in the laser-supported rocket the pressure has 
a substantial influence on the absorption and radiation processes which directly affect the 
available enthalpy of the propellant and the heat loading on the thruster wall. Figure 21 
shows the isotherms for a 300 kPa inlet pressure thruster, rocket No 4. The absorption 
coefficient and optically thin thermal radiation for a hydrogen plasma are approximately 
proportional to the 1.5 power of pressure. With a higher absorption coefficient in this 
plasma, compared to those of rocket Nos. 1-3, only a very small amount of the laser power 
directly escapes from the thruster (see Table 2); however, the optically thin radiation from 
the plasma also increases. This results in a low conversion efficiency from laser power to 
available thermal power of the propellant. The higher absorption coefficient also results in 
a plasma which is positioned further upstream against the laser beam. Two local tempera- 
ture maxima were predicted for this thruster. This results from two competitive processes 
which determine the local laser intensity: the attenuation by plasma absorption and the 
increase due to the converging laser beam. References 7-8 contain a detailed discussion of 
this subject. 


The vacuum specific impulse variations with area expansion ratio for the four calcu- 
lated thrusters are shown in Figure 22. The most advanced chemical rocket has an Isp 
of about 500 secs, and the calculated Isp for the four laser supported rockets are between 
700 and 1020 secs. At the exit of these rockets the ratio between thermal energy and 
the kinetic energy of the propellant is still relatively high, and further expansion of these 
rockets will greatly enhance the rocket performance. Although the propellant within the 
expansion section is expected to be far away from chemical equilibrium, the reasonable 
Isp is expected to be around 1500 secs for higher expansion ratios using rockets Nos. 1 
and 3. The highest Isp among the test cases is 1022 secs from rocket No. 3. As discussed 
before, this thruster has the propellant almost fully heated at the throat. Thruster No. 1 
has a little larger throat than rocket No. 3. and has a very small cold propellant stream 
near the throat. The calculated Isp of this rocket is only 20 secs less than thruster No. 3. 
With further increase of throat size, the size of cold stream propellant increases, and the 
calculated Isp decreases rapidly. The Isp for thruster No. 2 is an example. Rocket No. 4 
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Figure 21. Isotherms for rocket No. h (starting from 2,000 K to 15,000 K 
with 1,000 K increments) 
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AREA RATIO 

Figure 22. Vacuum specific impulse vs. area expansion ratio 



operated at 300 kPa and the plasma occurs further away from the throat than the lower 
pressure thrusters. This results in a larger cold propellant stream around the throat and 
in a much lower Isp than the lower pressure thruster. 

The heat loading on the inner wall of a good laser thruster should come from the 
plasma radiation, and a cold stream of propellant neax the inner surface will reduce the 
heat flux through the wall. This heat flux is used to heat the propellant from its storage 
enthalpy up to the inlet enthalpy. The irradiation along the wall of the calculated rockets 
are shown in Figure 23. The irradiation first increases from the end of the thruster to where 
the hot plasma is located, and then decreases monotonically along the axial distance for 
all test rockets. The maximum irradiation on the wall for all the rockets is almost the 
same, although the 300 kPa rocket has much higher total irradiation flux on the wall (see 
Table 2). In order to estimate the radiation heat loading on the wall, the emissivity and 
absorptivity of the wall surface are needed. If absorptivity is unity, the irradiation will 
be equal to the wall heat loading. Generally, the absorptivity is less then one, and the 
radiation heat loading will be smaller and spatially smoothed. The calculated maximum 
irradiation represents am upper limit of the radiation heat loading, and all maxima are 
much less then the heat loading for chemical rockets. 

Heat loading on the wall is also related to the thruster inlet propellant temperature 
which is regeneratively heated from the storage temperature. The propellant inlet temper- 
ature will influence the rocket performance as long as the cold stream of propellant exists 
near the wall at the nozzle. Rockets Nos. 1 and 3 have a very limited amount of unheated 
propellant at the throat, and the Isp for these rockets should not be strongly influenced 
by the inlet temperature. On the other hand, higher or lower inlet temperature will either 
increase or decrease the expected Isp significantly for rocket Nos. 2 and 4. An inlet tem- 
perature of 1500 K was used for the test rockets, and the expected inlet temperatures vs. 
three recovery efficiency of irradiation are tabulated in Table 3. Low pressure rockets Nos. 
1 and 2 need the recovery efficiency of 30 and 52 percent, respectively, for the adopted 
inlet temperature, 1500 K. If the irradiation is fully converted to the propellant enthalpy 
through regenerative heating for these two thrusters, the inlet propellant temperatures 
will be 2810 K and 2390 K, respectively. For rockets Nos. 3 and 4, the inlet temperatures 
will be around 3400 K, which is not realistic in a practical device, and probably a radiator 
would be needed to dissipate the wall heat flux and lower the inlet propellant temperatures. 

VH. EFFECTS OF CHEMICAL KINETICS 

The developed computer code for the laser-supported thruster calculation only consid- 
ers the chemically equilibrium conditions. The adopted numerical method is not computar 
tionally efficient for highly expanded supersonic flows either. The effect of finite rate chem- 
istry on these rocket performance for maximum area expansion ratio of 25 was evaluated 
using the TDK code 20 . The adopted kinetics equations axe the the same as that considered 
by McCay et al 21 for one-dimensional analysis. Six different options of the TDK code were 
considered in this study including two-dimensional equilibrium expansion, two-dimensional 
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Figure 24. Effect of finite rate chemistry on rocket performance 


frozen expansion, one-dimensional equilibrium expansion, one-dimensional frozen expan- 
sion, and one-dimensional kinetic expansion using two sets of rates compiled in Ref. 21. 
We experienced several difficulties during the practice of TDK calculations. For example, 
the TDK option of the TDK code fail to achieve a converge solution. It is because the fast 
chemical reaction rate under plasma conditions produces the numerical stiffness. 

In order to extract the output from the developed code of laser thruster calculations 
to fit into the required format for TDK code, some reasonable assumptions are needed. 
The predicted temperature and velocity profile across the throat (normally 61 grid points) 
from the developed code was divided into 11 stream tubes. The stagnation enthalpy for 
each stream tube was calculated. Then this information and chamber stagnation pressure 
were fed into TDK code with different selected options for rocket performance calculations. 

An example of this study for rocket no. 1 is shown in Figure 24. Since the predicted 
radial temperature profile at throat from the developed code is highly nonuniform and 
only eleven stream tube were used in the TDK code calculations, some errors may be 
introduced. The estimated inaccuracy for the above procedure in predicting the Isp value 
will be around 5 to 10 percent. The calculated values from two-dimensional options are 
always higher than the value calculated from one-dimensional analysis. The predictions 
considering finite rate chemistry do not differ too much from the equilibrium calculations. 
It is because that a large fraction of propellant at throat is in low temperature (between 
2,000K to 8,000 K) and is not influenced by finite rate chemistry. The Isp from this thruster 
is bounded by 1350 and 1550 secs from TDF and TDE calculations. 

Vm. SUMMARY AND CONCLUSIONS 

The computer code has been developed for the laser-sustained plasma flows. The 
adopted numerical algorithm has demonstrated the capability to handle complicated flows 
within a thruster. The physical models have been evaluated using existing experimental 
data. Through the parametric study of LSP in a constant area pipe, the power conversion 
efficiency from laser to the propellant and plasma positions were found can be controlled 
by adjusting the forced convection level and optical arrangements. From this experience, 
the computer code has been extended for the realistic thruster calculations. 

The performance of several realistic thrusters powered by a 30 kW laser has been 
studied. The results show that the laser supported rocket can have much higher vacuum 
specific impulse than the traditional chemical rocket. The Isp was in the range of 700 
to 1020 secs for four calculated thrusters having an area expansion ratio of four. The 
irradiation from the plasma to the thruster wall was also calculated, and the estimated 
maximum heat loading on the thruster is within the range of the chemical rockets. Among 
four calculated rockets, rocket No. 1 was the best practical rocket. This rocket has 
nearly the highest vacuum specific impulse among the test rockets, and both the estimated 
thruster heat loading and the inlet propellant temperature are reasonable. The estimated 
Isp for this rocket at expansion ratio up to 25 were also considered, and the bound of Isp 
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is from 1350 to 1550 secs considering the frozen and equilibrium expansion, respectively. 
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Numerical Method 




APPENDIX A 


A finite volume based computational scheme for a body-fitted coordinate sys- 
tem has been applied to the steady-state, axisymmetric transport equations in con- 
servative form. The method employs a nonstaggered grid such that the velocity 
components and all variables are defined at common grid positions. A general 
transformation was used to transform the governing equations from the nonorthog- 
onal, body-fitted, physical grid system to a new orthogonal grid mesh. 

The development of the numerical scheme closely follows Rhie and Chowt 13 ’ 14 !. 
A generalized differential equation in cylindrical coordinates will be transformed to 
the computational grid and the results then applied to the specific conservation 
equations. 

GENERAL FORMULATION 

The generalized transport equation may be written as 

^(p4>) + div(pV <f) = div(T 4 , grad <f>) + R<j>{r,x) (A.l) 

For steady-state conditions and axisymmetric, cylindrical coordinates, equation 
(A.l) becomes 


|^ u *> + + + rR *( r < x ) W 


A general transformation to the orthogonal plane may be employed by intro- 
ducing the new variables £ and rj (see Fig. 2) such that £ = £(r, 2 ) and rj = rj{r,x). 
It follows that the partial derivatives of any function / axe transformed by the 
following relations. 


fz = (r„/* - r*/„)/J 
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where J is the transformation Jacobian. 



J = x*r„ - x„r* 

(A.5) 

Introducing the relationships, 


U = ur„- VXr, 

(A.6) 

V = vx$ — vr$ 

(A.7) 
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ffi = x n + r n 
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92 = x^Xr, + r$r n 
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(A.10) 

produces the generalized equation in £ and tj coordinates. 


—(rpU*) + ±£(r,V,) = ~ - ft*,)’ 



(A.11) 


Integrating the transformed equation over a control volume shown in Fig. 2(b) gives 
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where f represents the direction of rotation about the axis of symmetry. Carrying 
out the integration yields 
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Let A rf = AC = 1 in the transformed plane for convenience and multiply through 
by J to get the finite difference approximation form for the integral conservation 
equation. 

[rpU4>Y w + (rpV^)" = - q 2 <f>ri ) 
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In the preceding equation, the small letters e,w,n,s represent the control vol- 
ume interfaces. The flow properties there must be interpolated from the known 
values at surrounding nodes denoted by capital letters. Assuming a piecewise-linear 
profile for <j> between the nodes yields the following approximations. 


(rpU</>) e « (rpU) e -(<j>E + <t>p) 

(A. 15) 


(rpU <f>) w « + 4>w) 

(A. 16) 


(rpV <j>) n « {rpV)n^[<j> N + <f> P ) 

(A.17) 


(rpV^,«(rpV),i(^ P + * s ) 

(A. 18) 


Also, 
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(^g ~ <£f) 
e A £ 


(^F ~ <f>w) 

e *£ 


(#i V - #f) 

A>? 


(#F ~ <fts) 

. A* 


(A.20) 


(A. 21) 


{A. 22) 


Values for [rpTJ), (rpV), (T^rqi/J) and, {T^rq^/J) are obtained using a hybrid 
scheme which will be discussed, in the next section. The value of 1/2 which arises 
in equations (A. 15) - (A.18) is the result of locating the cell interface midway 
between nodes. 

Making the appropriate substitutions into equation (A.14) yields 


i~Y~) M e + ~ (~Y~)w(4 > P + 4>w ) + C~Y~) n [<t>N + <f>P ) - (“““ ) a {<t>P + <f>s) 


- {4>E-4>p)~ ( < f > P-<t>w)+ {<f>N~<f>p) - ^ - -^ g3 ^ {<t>p-<t>s) 

(E^^ e )”] +rJSt(r,,t) (A.23) 

By taking all terms involving <f>p to the left hand side of the equation and all 
other terms to the right hand side, a relationship between <f>p and its neighbors is 
obtained. 


Ap<l>p = Ae4>e + Aw<t>w + -A-n^n + A.s^>s 
_[(L02^)^+(£ip^)"] +rJS^ v ,() [AM) 


A-4 



where 



(A.25) 


(A.26) 


(A. 27) 


(>1.28) 


Ap = Ae + Aw + An + As + ( rpU) e — (rpU) w + ( rpV) n — {rpV) a (>1.29) 


The quantities contained within the bracket in equation (A.24) result from the 
cross derivatives in the diffusion terms for the nonorthogonal coordinate system. 
Since these terms are small, they will be treated as known quantities and included 
with the source term. They may be evaluated through the following approximations. 

w ("t/ 92 ) 4&NB ~ <t>SE + <f>N - <f>s) (A.30) 

W ^ — ~ 4>sw + <t>N - <f>s ) (A.31) 

(^~y - <f>^j « (~J 4 {&NE ~ 4>nw +<f>E~ 4>w) (A.32) 

w ( t/ 92 ) (A.33) 

Again, terms such as (r ^rq^/ J) are evaluated by linear interpolation between nodes. 



In transforming the source term Ify(r, x) to the orthogonal coordinates 5 ^( 17 , £) , 
equations (A.3) and (A.4) should be applied in order to obtain the correct form. The 
source term may also be split to account for a linear dependence on the dependent 
variable <j>. 


5 = Sq + S P <t> p P (A.34) 

In equation (A.34), S is the average value of the source term through the control 
volume, Sc is considered to be the constant part of S, and Sp becomes a coefficient 
of <f>. Equation (A.24) takes the following form. 


Ap<f>p = Ap<j>E + + -A-n^n + A.s<f>s + Sc + Sp<f>p (A.35) 

Here, Sc now contains the terms involving the cross derivatives in the diffusion 
terms for the nonorthogonal coordinate system in addition to the constant part of 
5, and Sp must be negative to insure stability and to prevent physically unrealistic 
solutions. 

HYBRID SCHEME 

Special attention is also required concerning the evaluation of the A neighbor 
coefficients. Central differencing may be grossly inaccurate for these terms due to 
the non-linearities in the transport equations and because large convection terms 
mean that the ^-profile is far from linear between between nodes. 

A suitable formulation for evaluating these coefficients is the hybrid scheme 
which is identically the central difference formulation when the cell Reynolds number 
is less than 2 but reproduces the first order upwind scheme for a cell Reynolds 


number greater than 2. See Patankar^ 16 ! for a detailed discussion, 
scheme produces the following formulations. 

The hybrid 

“[( r 7).,('f)J-( r f). 

(A. 36) 


(A.37) 


(A.38) 
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(A.39) 


STABILITY 

Stability considerations for the discretized equation must now be examined. 
From equation (A.29), it is apparent that the Ap coefficient depends on its neighbor 
coefficients and the net mass source over the control volume. At steady state, the 
net mass source should be zero in a converged solution, but during the iteration 
process this may not be true. Also, a situation could arise such that all the neighbor 
coefficients are zero. Thus, the combination of terms which determine A p could at 
times become zero causing the finite difference solution for <j> to become singular. 
Clearly, this must not be allowed to occur. 

From inspection of equation (A.29), it is evident that, if the net mass source 
is positive, mass is being lost from the control volume and is carrying the property 
<f> out of the cell. In this case, the net mass source may be set to zero, and a new 
term added corresponding to the mass flow rate out, carrying the old property (<^p). 
For a negative net mass source, no adjustment is necessary. Introducing the above 
modification yields a new discretized equation. 

Ap<f>p = Ae4>b + Aw<t>w + An<I>n + As<f>s + Sc + Sp<f>p (A.40) 


Ap = Ae + An + Aw + As (A.41) 


S c = max {(rptf)e - (rpU) w + (rp7) n - (rpV).},0 <f>* P + S% 9 (A.42) 


S P = -max^{rpU) t - ( rpU) w + ( rpV) n - (rpY).J,oJ + S° p rg (A.43) 

The terms S^ 9 and S^ a correspond to the original source terms in equation (A.35), 
and 4>* p corresponds to the previous iteration result for the dependent variable. 
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RELAXATION 


The above algorithm will guarantee the convergence of linear differential equa- 
tions, but this criterion may not be applied for nonlinear differential equations. To 
prevent divergence and to slow down the changes between iterations, underrelax- 
ation of the dependent variable is employed. The discretization equation (A.40) is 
first reformulated to give 


[Ap — Sp)<f>p = ^ A<f> + Sc 


(A.44) 


where the sum of the product of the neighbor coefficient and dependent variable 
has been represented by the summation symbol. Dividing through by (Ap — Sp), 


<t>p = 


13 A0 + Sc_ 
Ap — Sp 


(A.45) 


Now, if <f>* p is denoted as the old result from the previous iteration and added to 
and subtracted from the RHS of equation (A.45) the result is 

= (^71 + S^ “ ( A ' 46 > 

Introducing a relaxation factor a gives 

<f>P = <t>* P + a “ <*p) (A.4 7 ) 


Multiply through by (Ap — Sp). 


(Ap — Sp)<f>p = (Ap — Sp)<f>p + a 


yi A<f> + Sc — (Ap — Sp)<j>p 


(A.48) 


Finally, divide through by a and rearrange to get 


{Ap ~ Sf) H + 5c + (1 - «) (Ap ~ Sf) r P (AA9) 
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This is the new discretization equation to be solved for <j>p. The general formulation 
has been completed and the results may now be applied to the governing transport 
equations. 

EQUATIONS OF MOTION 

The steady-state, laminar Navier-Stokes equations and the steady-state conti- 
nuity equation are solved to determine the fluid motion. The momentum equation 
is expressed in general transport form and solved for the velocity components. The 
continuity equation, however, must be reformulated into a pressure correction equa- 
tion which becomes appropriately coupled to the momentum equations. 

The solution procedure used to solve these pressure coupled equations is a 
marriage between the algorithm developed by Rhie t 14 l and the PISO method of 
Issal 16-18 !. The basic idea is to correct the mass flux obtained from the solution 
of the momentum equations such that the discretized pressure equation becomes a 
transport equation for compressible flow instead of a Poisson equation. The split 
operator concept of PISO is utilized to solve the pressure correction equation. 

The development of the discretized momentum equations begins by writing 
them in cylindrical coordinates. 


d . . d . . d . du . d . du. dP , . 

^(r^u) + -(r/nm) = ^(rn—) + ^(r^) - r— + rfi«(r, i) 


(A.50) 


d , \ , d . . d . dv. d . dv. dP _ . . 

^(rpuv) + -(rpvv) = + &&■§;) ~ + r. R,(r,x) 


(A.51) 


where the terms involving the pressure gradients are a part of the source term. From 
comparison with equation (A.2), the parameters for the ^-direction momentum are 
<f> = u and T = n, and the parameters for the ^-direction momentum are, in a 
similar manner, <f> = v and T = fi. 

Transforming to the orthogonal grid yields. 


Apup — ^ Au 


4. ( ^L3l du \ n 
A J dr,) 9 + \ J dZ),. 


+ 


dr_dP_ dr_d£ 
T dr} d£ T dZ drj 


+ rJ5„(r7,$) 


(A.52) 
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- V- , r frnavV . 

Apvp= ^ Av -[{—^). + \—w.. 

-j + rJS v (r), £) 


dx_dP__ dx dP 
T drj d£ T d£ dr) 


(A.53) 


where the terms inside the brackets axe considered constant and taken as part of 
the source term. 

Following the formulation of equation (A.40) and employing the underrelax- 
ation technique previously discussed, equations (A.52) and (A.53) take the following 
discretized form for the dependent variables u and v respectively. 


(Ap ~ Sp) u F = X> + Sc + (1- a) (Af ~ Sp) u- P (A. 54) 

(Ap-g P ) ^ = ^ A „ + 5c + (1 _ tt) (Ap-Sp) p . {AM) 

Note that the coefficients and source terms in the above relationships are unique 
with respect to the dependent variable. The underrelaxation factors need not co- 
incide either. The dependent variables u and v axe solved for from the above 
discretization equations. 

PRESSURE CORRECTION EQUATION 

As previously stated, the split operator concept of Issa I 16-18 1 J as developed 
in the algorithm of Rhie l 14 ’ 15 ), is implemented in this work to solve the coupled 
momentum and continuity equations. The fundamental approach is to correct the 
preliminary mass flux through the following perturbation relationship. 


rpGi = r(p * + />')(<?< + * rp'G \ + rp'{G\ + Gj) (A.56) 

The variable G, denotes the contravariant velocity, and the superscripts *, ' rep- 
resent the intermediate iteration values and the correction increment, respectively. 
The perturbation relationship may then be combined with the continuity equation 
through a series of split operations producing a sequence of pressure correction 
equations in transport form; thus, the correction equations may be discretized in 


i 
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the same manner as the previously developed method for the general, transport 
differential equation. The solution of these correction equations gives velocity and 
pressure fields which satisfy both the momentum and continuity equations. An 
important characteristic of this method is the ability to handle both incompressible 
and compressible flow regimes over the entire flow field. A formal derivation of the 
procedure follows. 

First, recall relations (A.6) and (A.7) which relate the contravariant velocity 
components to the physical velocity components and substitute for u and v from 
the momentum equations (A. 52) and (A.53). The results are 


u = 

|fT € (u,t;) - B^A^P - C^Ar,P + St 

(A. 57) 

. . dr Y A.U dx Y Av 

H ^=a, A f ~ 3n 

(A.58) 


_ r dr dr r dxdx 

* A “ drj dri + A” dr} dr\ 

(A.59) 


_ r dr dr r dxdx 

* A“ d£ dr} A" d£ drj 

(A.60) 


= dr Su dx S 

dr} A* drj A« 

(A.61) 

V = 

v ) “ — CrjA^P + Sfi 

(A.62) 

tt (u v \ _ dx Y,Av drY.Au 
~ Bi A’ 9( A 9 

(A.63) 

_ , . r dx dx r dr dr 

~ ~A” p ~dl~dl + ~A*~dl'dl 

(A.64) 

~ . r dxdx r dr dr 

” J A* p dr,dt A^drjdt 

(A.65) 
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(A.66) 


B dx S v dr S u 

n ~d(A°~diA} 

The A i symbol represents the finite difference equivalent of d/dx{. For ease of 
presentation, equations (A.57) and (A.62) may be written in indicial notation as 


Gi = 


Hi(u, v) - BiAiP - CiAjP + Si 


(A.67) 


If a pre liminar y pressure field is imposed, the mass fluxes across the control 
volume boundary are established yielding an initial predictor equation. 


r P*G\ = rp* 




D a n» 


(A.68) 


Equation (A.68) is effectively a momentum equation for the contravariant velocities 
U and V . The P* pressure field, however, will not normally allow the correspond- 
ing continuity equation to be satisfied throughout the flow field. The essence of the 
split-operator concept is to establish momentum and continuity satisfying velocity 
and pressure fields through a sequence of split operations using the perturbation re- 
lationship. Instead of correcting the pressure with one operation as in SIMPLE, new 
velocity and pressure fields are evolved through a succession of correction operations 
which attempt to update the field values in a step-wise manner. The procedure is 
to apply the perturbation relationship (A.56) - neglecting the product of the two 
perturbation variables - such that three correcting operations are formulated as 
follows: 

(a) First Correcting Operation 


rp**G* { * » r(p ** - p*)G\ 


+rp* 


Hi{ « V) - BiAiP** - CiAjP* + Si 


(A. 69) 
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(b) Second Correcting Operation 


rp G f » r(p - p JGi 


+rp* 


Hi{u\v *) - BiAiP*** - CiAjP** + 5, 


(A.70) 


(c) Third Correcting Operation 


rp****G*i*** a r{p**** -p*) 


+rp* 


v***) - B^P **** - CiAjP ** + Si 


(A71) 


The number of * superscripts denotes the intermediate iteration sequence. Note 
that several alternative approaches exist for operator-splitting and that the method 
presented here is only representative of Rhie’s ( 14 1 algorithm. Another important 
point to remember is that each correcting operation acts as the predictor for the 
proceeding operation. There are three unknowns in each of the above equations - 
that is, the corrected velocity, pressure, and density; fortunately, the equation of 
state allows the density to be expressed as a function of pressure and temperature. 
If the temperature is assumed constant over each complete iterative cycle and then 
updated at the end of each cycle, the correcting operation equations contain only 
the corrected velocity and pressure as unknowns. Supposing that the new pressure 
field is obtained through some other means indicates that the corrected velocities 
may be calculated, explicitly. This is exactly the approach to be taken here. 

The formulation of the desired pressure correction equations is accomplished by 
subtracting the predictor equation from its respective correcting operation equation, 
applying the equation of state, and enforcing continuity for the unknown corrected 
velocity. The density increment p\ as expressed as a function of P' by the ideal gas 
law, is underrelaxed as follows: 



(A.72) 


/? is the underrelaxation factor; R is the gas constant; and T* is the static tem- 
perature from the previous iterative cycle. The formulation procedure yields the 
pressure correction equations which follow: 
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(a) First Pressure Correction 


A.^G;(P" - />•) - A ir p'B t A t (P" - P-) = -Atrp-G- (A. 73) 


(b) Second Pressure Correction 


Aij^G-(P— - P") - Asp-BiAtiP— - P") 


= Airp‘CiAi(P" - P‘) - Airp"Gy 


{A. 74) 


(c) Third Pressure Correction 


Ai RF G *( P **** ~ " A ‘ r / J * 5 * A *( P **** “ -P***) 


= -A t V 


.&(«***,»;***) -H{U\v*) 


- A irp***Gt** 


(yi.75) 


These correction equations are expressed in the form of the general differential 
equation. The first team on the LHS represents the convective part; the second term 
on the LHS represents the diffusive part; and the entire RHS is the source part. The 
correction equations are effectively the transport equations for compressible flow. 

The physical implication of these equations is obtained by considering the grid 
cell Reynolds number. For convenience, let’s assume that an orthogonal grid cell 
exists with a u component of velocity only. The grid cell Reynolds number for such 
a situation would be 


Re e = 


convection 


= Ptzz = = uM2 


diffusion RT 


Employing an upwind-difference scheme implies that the critical grid cell Reynolds 
number may be made to coincide with M = 0.3 by adjusting the value of w. Such a 

method allows the convective terms to vanish in the incompressible flow limit and 
the same terms to dominate when the Mach number increases and the compressibil- 
ity effects become consequential. For a nonorthogonal grid and multiple components 
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of velocity, the procedure is more complex, but the principle and physical correspon- 
dence still apply. At increasingly higher Mach numbers, artificial viscosity must be 
introduced to increase the smoothing properties of the correction equations and to 
suppress numerical oscillations. The term used by Rhie ( 14 1 as taken from Shamroth 
et al t 21 l may be implemented in the present work in the form 


Hart = max 



p£ 



(A. 77) 


Because a nonstaggered grid is used in this work, an oscillatory pressure field 
could result. To eliminate this possibility, the G\ value at the cell boundaries 
when evaluated for the RHS of equation (A. 73) includes a fourth order pressure 
dissipation term. These flux terms are evaluated according to Rhie t 14 l in the form 


u; +m = (v; + v‘ i+l )i2.o+ Ci (p; +m - s p; n + 3 p; - (x.rs) 

The coefficient Ed is a function of Ap and the metric coefficients. The dissipation 
term allows the linearized flux and pressure gradient to become coupled. 

The pressure correction equations are discretized in a form similar to the gen- 
eral transport equation. The pressure correction increment is solved for and sub- 
stituted back into its respective correcting operation equation to determine the 
corresponding velocity components. Once this is done for all three correction equa- 
tions, the energy transport equation is solved and all properties updated. The entire 
sequence is then repeated to convergence. 

ENERGY EQUATION 

The energy transport equation requires special consideration since a high tem- 
perature plasma is to be modeled. Radiation absorption from the laser power source 
and radiative heat loss must be incorporated into the physical model as described 
in the main body of the report. 

Assuming that the flow is laminar, steady, axisymmetric, and has variable 
physical properties, the energy transport equation may be written in cylindrical 
coordinates as 


A-15 


(A.79) 


d . d . . d 

—(rpuh) + —(rpvh) = - 


(K,n dh 
\ c p dx 


) 


+ 


d 

dr 



aIi-q ra d + rR{r,x) 


The last three terms are considered known from the previous iteration and comprise 
the total source term. 

Prom comparison with the general transport equation, the dependent variable 
is seen to be <f> = h and the diffusion coefficient is T = K e ff/C p where K e ff/C p is 
compiled in a table as a function of pressure and temperature and is reevaluated at 
each iteration. 

If equation (A.79) is transformed to the orthogonal grid, discretized, and formu- 
lated in the same manner as equation ( A.40) , the following form for the dependent 
variable h results. 


{Ap ~ Sp) h P = Y,Ah + Sc + ( 1 - «) Sp) h' P (A. 80) 

BOUNDARY CONDITIONS 

In the present numerical code, all boundaries coincide with the control volume 
sides. For the Dirichlet condition, the prescribed value of the dependent variable 
is specified at the control surface boundary; the correct convection and diffusion 
coefficients are calculated, and a new total coefficient and source term are obtained. 
For the Neumann conditions, the normal derivative to the boundary is specified 
by setting the appropriate neighbor coefficients to zero and marking the prescribed 
boundary flux into the source term. The specific boundary conditions for the present 
application axe presented in the main body of the papers. 

SOLUTION PROCEDURE 

Due to the immense size of the matrices representing the systems of equations, 
direct solution techniques were considered impractical, and an iterative method 
was employed. Because a node point by node point Gauss-Siedel iterative scheme 
has slow convergency rate, a line by line method using the TVi-Diagonal Matrix 
Algorithm (TDMA) was implemented. A complete description of this method is 
given by Patakar t 15 I. It is important to remember that the iterative method requires 
several passes through the flow field to approach convergence. Of course, the system 
need not converge completed for intermediate steps, but it should converge enough 
to provide sufficiently reasonable value for the next step toward the final solution. 
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Abstract — Laser thermal propulsion refers to a concept in which power from a laser is beamed to a vehicle, 
where it is absorbed and used to heat a propellent gas that produces thrust at high values of specific 
impulse. The success of the concept depends upon the ability to create laser sustained plasmas that are 
stable in the presence of a forced convection flow and that absorb essentially all of the incident laser 
radiation. In addition, the thermal radiation losses from the plasma must be controlled to realize high 
performance without the penalty of thermal radiators. Extensive experiments have been performed and 
analyzed utilizing plasmas sustained in flowing argon by a carbon dioxide laser of less than 1 kW incident 
power. Spatially detailed measurements of the temperature in the plasmas have resulted in a detailed 
understanding of the energy conversion processes occurring within the plasma. This understanding has 
provided an explanation of the limited range of pressure and power observed for stable operation of laser 
sustained plasmas and suggests the means by which that range can be extended. A computational model 
based on a Navier-Stokes description of the flow has been developed which utilizes ray tracing for the 
optical interaction and temperature dependent transport and thermodynamic properties for the plasma. 
This model has been verified using the detailed experimental results for the argon plasmas and has been 
extended to provide predictions for higher laser powers in hydrogen. 


I. INTRODUCTION 

Laser thermal propulsion refers to a concept in which 
power from a laser is beamed to a vehicle, where it 
is absorbed and used to heat a propellant gas that 
produces thrust at high values of specific impulse. 
This concept has the advantages that pure hydrogen 
can be used as the propellant and that the weight of 
the propulsion system can be small, since the power 
used to heat the propellent does not have to be 
carried aboard the vehicle. Several methods can be 
used to absorb the laser energy, but we shall restrict 
our discussion to the case where the energy is ab- 
sorbed through the process of inverse bremsstrahlung 
in a high temperature plasma of the propellent gas. 
When pure hydrogen is used as the propellent gas, 
specific impulse values in the range of 1000- 1500 s 
can be expected and significant improvements are 
obtained for the fractional payload that can be 
delivered from low earth orbit (LEO) to geosyn- 
chronous earth orbit (GEO)[l]. 

A schematic representation of the motor is shown 
in Fig. 1. The incoming laser beam is focused into the 
absorption chamber where it sustains a plasma at 
temperatures of approximately 16,000-20,000 K. In 
contrast to a normal combustion chamber, the high 
temperature plasma occupies only a small fraction of 
the chamber volume, and mixes with the incoming 
propellent to produce a flow having a bulk tem- 
perature of approximately 3000-5000 K that can be 

tPaper IAF-86-17S presented at the 37th Congress of the 
International Astronautical Federation , Innsbruck, Austria, 
4-11, October 19S6. 
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expanded through the nozzle to produce thrust. Due 
to the high temperature in the laser sustained plasma 
(LSP), a substantial portion of the absorbed power 
will be radiated from the LSP region and absorbed by 
the chamber walls. In order to utilize this radiated 
power effectively, the chamber walls are cooled by the 
incoming propellent in a regenerative cycle. If the 
power radiated from the plasma can be limited to that 
required to raise the propellent from storage tem- 
perature to an acceptable inlet temperature, then all 
of the absorbed power can be utilized to heat the 
propellent and no external radiator will be required. 
Sufficient mixing of the inlet propellent and the high 
temperature region within the LSP must be insured to 
produce the bulk temperature required for the nozzle 
inlet. 

The fraction of absorbed laser power that is lost 
from the plasma as thermal radiation has a profound 
influence on the performance that can be achieved by 
the thruster. If the acceptable inlet temperature limits 
the enthalpy of the inlet flow to /*<,, then the final bulk 
enthalpy will be limited to h T = V a RAD> where 
is the fraction of the absorbed laser power that is 
radiated to the walls by thermal radiation. The bulk 
enthalpy, and therefore the impulse of the thruster, 
will be determined by the radiated fraction of the 
absorbed laser power that may be absorbed by 
regenerative cooling. It is clear that the success of this 
concept rests on the ability to sustain a stable plasma 
in a forced convective flow which can absorb substan- 
tially all of the incoming radiation and limits the 
fraction of the absorbed laser power that is lost from 
the plasma region through thermal radiation. 
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Fig. 1. Schematic representation of a rocket motor that operates by laser thermal propulsion. 


Experimental studies of laser sustained plasmas 
have been conducted at The University of Tennessee 
Space Institute (UTSI) since 1983 to investigate their 
suitability for laser thermal propulsion. The objective 
of this research has been to understand the basic 
physical mechanisms of the LSP and determine their 
influence on stability, laser power absorption, and 
thermal radiation loss. The approach employed was 
to create the LSP in flowing argon at moderate 
carbon dioxide laser powers (less than 1.0 kW) and to 
obtain detailed measurements of the plasma tem- 
perature field using a modem digital image pro- 
cessing computer to acquire spectral images of the 
plasma emission. Using the measured temperature 
field, the laser power absorption and thermal plasma 
emission could be calculated at any point within the 
LSP to provide a detailed understanding of the 
energy conversion processes[2]. 

A computational model for the LSP in a forced 
convective flow has also been developed by Jeng and 
Keefer(3] to aid in the interpretation of the experi- 
mental results[4] and to provide a predictive capabil- 
ity for LSP’s operating at higher power and in 
hydrogen[5]. This model for the LSP incorporates the 
temperature dependent thermodynamic and trans- 
port properties for the argon plasma within a Navier- 
Stokes model of the fluid mechanics. In addition, the 
model incorporates geometric ray-tracing for the 
optical configuration that permits the use of laser 
beam profiles of arbitrary power distribution and 
includes the effects of lens aberrations. The resulting 
model is solved numerically using a Spalding- 
Patanker SIMPLE algorithm. 

2. EXPERIMENTAL INVESTIGATION 

The primary objective of our experimental in- 
vestigations has been to understand the physical 
mechanisms that determine the characteristics of the 
LSP in the regime where forced convection dominates 
the flow. In most previous experiments the LSP has 
been operated in a large chamber[6], or in the open 
air[7] where the flow is dominated by free thermal 
convection. The stability of the LSP was of particular 
interest to us, since stability is of crucial importance 


for propulsion and because it has been a subject of 
some discussion in the literature[8]. Our approach has 
been to sustain plasmas in argon using a relatively 
low power (1000 W) C0 2 laser and to obtain detailed 
spectroscopic measurements of plasma temperature 
for a variety of combinations of pressure, flowrate, 
laser power and focusing geometry. The measured 
temperature fields were then used to determine the 
detailed absorption and radiation processes that oc- 
cur within the plasma. 


2.1. Experimental procedure 

The plasmas sustained in these experiments were 
contained in a 22 mm, inside diameter, cylindrical 
quartz flow tube which had a converging entrance 
section. A diagram of the experimental apparatus is 
shown in Fig. 2. 

The laser power was provided by an axial-flow, 
electric discharge, carbon dioxide laser which pro- 
vides a nominal 1 kW, continuous beam. The beam 
has an annular profile produced by an unstable 
oscillator of 1.5 magnification. The beam was ap- 
proximately 54 mm in diameter with good axial sym- 
metry, and was focused collineariy with the axis of 
symmetry of the flow tube. The apparatus was ori- 
ented with its symmetry axis in the vertical direction, 
with the laser beam propagating vertically upward, to 
avoid nonaxisymmetric buoyancy effects. This 
configuration generally produced a plasma with the 
good axial symmetry that is required for accurate 
diagnostics. 

The primary measurements acquired from the ex- 
periments were digital images of the plasma. These 
images were obtained by a CID video camera whose 
spectral range was limited by a narrow bandpass 
interference filter. The filter was chosen to isolate a 
1 nm wide region of the continuum radiation from 
the plasma, free of line radiation, and centered at 
626.5 nm. The digital images were acquired by a 
digital image processing computer, which allowed 
images to be acquired at video frame rates (60/s) and 
averaged into memory. The entire image acquisition 
system was calibrated through direct substitution of 
a standard lamp of spectral radiance. Additional 


I 


B-3 


Laser thermal propulsion using laser sustained plasmas 


369 




Exhaust 



Fig. 2. Schematic diagram of the experimental apparatus. 


details of the apparatus and experimental procedure 
have been described in earlier papers[2]. 

More than 125 different experiments have now 
been completed using lenses of two different focal 
lengths and a wide range of pressure, flow and power 
conditions. A new, transform domain method was 
developed to permit efficient and accurate Abel in- 
version of the large quantities of data that were 
obtained[9]. By applying appropriate calibration fac- 
tors, the emission coefficient was obtained as a func- 
tion of position in the plasma, and from this the 
temperature field within the plasma was obtained by 
interpolation in tabulated data published by 
Morris[10]. 

Detailed geometric ray tracing techniques were 
employed to calculate the spatial power absorption at 
each point in the plasma, using the measured tem- 
perature field to determine the laser absorption 
coefficient. To calculate power absorption within the 
plasma, the measured annular profile of the incoming 
laser beam was broken into discrete rays, each repre- 
senting a fraction of the incident radiation. Each ray 
was traced, using geometric optics, through the so- 
dium chloride focusing lens and window and then 
through a grid system superimposed on the measured 
temperature field of the plasma. The absorption 
coefficient at 10.6 /im was calculated according to the 
method described by Kemp and Lewis[U], using the 
thermodynamic data tabulated by Drellishak[12], and 
the Gaunt factors published by Karzas and 
Latter[13]. At each step through the grid, the power 


absorbed from each ray was subtracted from the ray 
and added to the plasma. This procedure gave both 
the spatial distribution of power absorption and the 
total power absorbed by the plasma. 

It was also possible, using the measured tem- 
perature fields, to calculate the power lost from the 
plasma due to the optically thin portion of the 
thermal radiation. The thermal radiation losses from 
the plasma were determined using the relation pub- 
lished by Kozlov for optically thin radiation[8]. 

2.2. Experimental results 

The objective of these experiments was to in- 
vestigate the effect of high incident flow rates on the 
behavior of laser-sustained plasmas. At flow rates of 
less than approximately 30 cm/s, the position and 
shape of the plasma relative to the incident laser 
beam was essentially unaffected by the small forced 
convective flow. The plasma stabilized in the region 
of the focused beam where the incident intensity was 
just sufficient to balance the thermal radiation, con- 
duction, and convection losses from the plasma. In 
general, this means that the plasma will stabilize at a 
position upstream of the focus that depends on the 
incident power, focusing geometry, and ambient pres- 
sure. The plasma tends to fill the beam at the position 
of stabilization, and therefore grows larger in dia- 
meter when it stabilizes at a position further up the 
beam. When the incident flow rate is increased 
significantly beyond 30 cm/s, convection becomes rel- 
atively more important in the energy balance and 
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Fig. 3. Isothermal contour plots of two plasmas sustained 
using the 203 mm (8-in.) lens and a chamber pressure of 
2.5 atm. The contour interval is 500 K, with the outer at 
10,500 K. 


tends to force the plasma to stabilize at a position 
further down the beam and closer to the focus. The 
plasma again tends to fill the beam, but because the 
beam is smaller in diameter near the focus, the 
plasma is also smaller in diameter and lies in a region 
of greater laser intensity. 

Figure 3 shows an isothermal contour plot for the 
high temperature region of two plasmas sustained 
under similar conditions of pressure and laser power, 
but with substantially differing incident flow speeds. 
The contour interval is 500 K with the outer isotherm 
at 10,500 K, and the focal point of the laser (defined 
as the location of the circle of least confusion) is at 
the location labeled 50 mm on the x-axis. The flow is 
from left to right, and the laser beam is incident from 
the left. The path of the incident laser beam through 
the plasma is also illustrated in Fig. 3 where the lines 
represent the approximate boundaries of the annular 
beam in an unrefracted path through the plasma. For 
the low-flow case, the major portion of the plasma is 
upstream of the laser focus, and there is a local 
off-axis temperature maximum. The off-axis tem- 
perature maximum in this case occurs near the an- 
nular focus of the beam produced by spherical aber- 
rations, as discussed in an earlier paper[14]. By 
comparison, the high-flow plasma stabilized with a 
significant fraction downstream of the focus, and has 
no significant off-axis temperature maximum. In ad- 
dition, the maximum temperature of the high-flow 
plasma is higher than that for the low-flow case; 
1 5,600 K compared to 14,600 K. 

The path length of the beam through the two 
plasmas shown in Fig. 3 is comparable, but in the 
high-flow case the path traverses, on the average, a 
region of higher temperature and, therefore, a larger 
absorption coefficient, particularly in the region of 
highest laser intensity near the focus. This results in 


a significantly higher fractional absorption of the 
laser poser; 78% compared to 67%. At the same 
time, because the overall size of the high-flow plasma 
is smaller than the low-flow plasma and because the 
highest temperature region is confined to a smaller 
volume near the axis, the overall thermal radiation 
loss decreased from 53 to 45% of the incident laser 
power. As a result of the different spatial character- 
istics of the plasmas, the conversion from laser power 
to flow enthalpy increased from 14 to 32% as a result 
of the increase in flow. 

The two sets of isotherms shown in Fig. 3 represent 
the endpoints of a series of different flow velocities. 
The data for fractional absorption of laser power and 
thermal radiation loss at the intermediate values are 
shown in Fig. 4. 

The absorption fraction passes through a max- 
imum near an incident flow velocity of 300 cm/s, 
while the fraction radiated reaches a maximum near 
100 cm/s. The difference between the absorbed and 
radiated fractions represents the fraction of the inci- 
dent laser power that is used to increase the enthalpy 
of the flow, and this quantity increases with flow up 
to a velocity of approximately 300 cm/s, after which 
it is nearly constant at 30%. 

The focusing geometry used to sustain the plasma 
also has a strong influence on the power absorption 
and the stability characteristics of the LSP. We have 
conducted a series of experiments using two different 
focal length lenses: 203 mm, corresponding to an 
f/number of 3.8, and 305 mm, corresponding to an 
f/number of 5.6. Isothermal contour plots of plasmas 
sustained with the two lenses at a pressure of 2.5 atm 
are shown in Fig. 5. 

Both plasmas were sustained at approximately 
720 W and 4 m/s incident flow velocity in argon at a 
pressure of 2.5 atm. The f/3.8 (203 mm) plasma is 
centered about 2 mm downstream of the focus, and 
has a peak temperature of 15,600 K, while the f/5.6 



Flow 

Fig. 4. Variation of the laser power absorption [▼] and 
thermal radiation loss [A] with incident flow velocity for a 
chamber pressure of 2.5 atm using the 203 mm focal length 
lens. The powers are shown as a fraction of the 720 W 
incident laser power. 
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Fig. 5. Isothermal contour plots for two plasmas sustained 
at 2.5 atm in a 4 m/s incident flow with 720 W incident laser 
power using 203 mm (8-in.) and 305 mm (12-in.) focal 
length lenses. The contour interval is 500 K and the outer 
contour has a value of 10,500 K. 


(305 mm) plasma is centered within 0.5 mm of the 
focal point, is just slightly smaller in diameter, and 
has a peak temperature of 16,200 K. Because of the 
narrower angle of convergence of the beam, the flow 
speed required to force the plasma into the down- 
stream region of the focal zone, where the efficiency 
begins to deteriorate, is greater for the larger 
f/number lens. This suggests that stable plasmas may 
be sustained using very long focal length lenses, 
provided that the forced convective flow is sufficiently 
fast. 

A comparison between the results of a series of 
experiments in argon at a pressure of 2.5 atm using 
both the 203 mm lens (f/3.8) and a similar series using 
the 305 mm lens (f/5.6), is shown in Fig. 6. Note that 
the fractional laser power absorption is greater for 



Flow (cm/s) 

Fig. 6. The effect of incident flow velocity on laser power 
absorption [V] and thermal radiation [A] for plasmas 
sustained at a pressure of 2.5 atm using 203 mm (8-in.) and 
30omm (12-in.) focal length lenses. 


the longer focal length lens at the higher incident flow 
velocities, but that it is smaller for the lower incident 
flows. It is also seen that the conversion of laser 
power to flow enthalpy is greater for the longer focal 
length lens at the higher flow velocities. For the 
longer focal length lens, the fractional laser absorp- 
tion and the fractional thermal radiation losses are 
clearly tending toward the same value as the flow 
decreases, indicating that the plasma will become 
unstable in this flow region. In an earlier paper[14], 
we reported that we were unable to initiate or sustain 
a plasma pressures greater than two atmospheres 
using the 305 mm focal length lens. The incident flow 
speed in that set of experiments was limited to less 
than 15 cm/s. As the pressure increased, at a constant 
mass flow rate, the plasma would move upstream 
until, at 2 atm, it became unstable. At the much 
higher flow speeds used in the later experiments, we 
had no difficulty initiating and sustaining plasmas 
with the 305 mm lens at pressures up to 3 atm. These 
results suggest that earlier experimental maps for 
regions of stability of the LSP[6] are not unique, but 
depend on the particular optical geometry employed. 

The plasma position and shape also depend 
significantly on the pressure and incident laser power. 
At higher chamber pressures, the absorption 
coefficient is sufficiently large to allow the plasma to 
move upstream into a region of lower laser intensity. 
With the annular laser beam this can, in extreme 
cases, result in large diameter plasmas with strong 
off-axis temperature maxima. Because of the lower 
laser intensity in the annular focus relative to the 
axial focus, the magnitude of a typical off-axis tem- 
perature maximum is lower than one on the axis, and 
this results in a lower fractional power absorption. At 
the same time, the large plasma diameter ensures that 
the thermal radiation losses remain large. If the 
incident flow velocity is insufficient to hold the 
plasma near the focus, then the plasma will attempt 
to move far up the beam away from focus and 
become unstable. 

The effect of incident laser power on the position 
and shape of the plasma is similar to that of pressure. 
For a given incident flow velocity, an increase in laser 
power causes the plasma to move up the beam to a 
new position of stability and a larger plasma di- 
ameter. If the plasma moves too far up the beam, the 
thermal radiation loss increases and, at the limit of 
stability, the radiation losses exceed the absorbed 
laser power. These positional effects would explain 
both the power instability and the pressure instability 
observed by Generalov et ai [6]. 

Detailed analysis of a large number of laser sus- 
tained plasmas has revealed the complex interactions 
between the flow and the energy conversion processes 
within the plasmas. These experiments also suggest 
that stable laser sustained plasmas may be operated 
over a wider range of conditions than previously 
reported by using appropriate combinations of opti- 
cal geometry, pressure and incident flow velocity. 
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3. THEORETICAL ANALYSIS 


3.1. Numerical model 


The first theoretical model for the laser sustained 
plasma was developed by Raizer, who proposed a 
one-dimensional model based on an analogy to a 
subsonic combustion wave[15]. This model was ex- 
tended by Jackson and Nielsen to include radiative 
transport 16], and by Kemp and Root to calculate 
hydrogen plasmas for laser propulsion[17]. Keefer et 
al. applied a different numerical method for solving 
the one-dimensional model that eliminated a problem 
with numerical convergence[18]. A quasi-two- 
dimensional model was proposed by Batteh and 
Keefer in which the flow and the laser beam were 
both assumed parallel to the axis and simple assump- 
tions were made for the plasma properties to permit 
a closed form solution[19]. Glumb and Krier ex- 
tended the quasi-two-dimensional model to include 
focusing of the laser beam and used a numerical 
solution to permit the inclusion of realistic thermo- 
dynamic and transport properties, but retained the 
axial flow, constant mass flux assumption[20]. 

A new model has been developed by Jeng and 
Keefer that is based on a laminar Navier-Stokes 
description of the flow and which utilizes geometric 
ray tracing to describe the optical interactions in the 
plasma[3]. This model has been verified by com- 
parison with the detailed experiments described 
above[4], and has been extended to predict the prop- 
erties of hydrogen at higher laser powers[5]. 
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Fig. 7. Details of the solution domatr ujeu tor the numerical 
calculations. 



Ideally, the model would simulate the flow geome- 
try shown in Fig. 1. However, in order to gain more 
understanding of the interactions between laser and 
plasma, the simple flow geometry as shown in Fig. 7 
was adopted. It is assumed that this domain is 
axisymmetric, and a cylindrical coordinate system 
was employed. The incoming working fluid enters 
from the bottom, and the originally collimated 
10.6 pm carbon dioxide laser beam is focused into the 
gas by an appropriate lens. 

It was assumed that the flow is laminar, steady- 
state, axisymmetric and has variable physical proper- 
ties. The pressure of the flow is relatively high, so that 
local thermodynamic equilibrium (LTE) can be as- 
sumed with little error. Therefore, the plasma can be 
described by a single temperature, and its physical 
properties are only a function of this temperature and 
pressure. Thermal radiation was divided into two 
parts: an optically thin portion in which all radiation 
escapes from the plasma; and an optically thick 
portion which can be described as a diffusion process. 
Since the mach number for the flow is small, the 
kinetic energy and viscous dissipation were neglected 
in the energy equation. Following the above assump- 
tions and using incompressible flow viscous terms, 
the equations of conservation of mass, momentum 
and energy for the flow can be written as: 


d(pu) 1 d(rpv) 
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where /, is the local laser intensity described by the 
following equation 


dli 

dr 


-a/,, 


(5) 


where s is the distance along the laser propagation 
direction. The laser beam path through the optical 
system was calculated by geometric ray tracing. Re- 
fraction of the individual rays due to the in- 
homogeneous refractive index within the plasma were 
neglected in the calculation. Since Welle et al\ 2] had 
found that this assumption produced little error in 
the experimental results. Diffraction of the beam 
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when combined with lens aberrations can have a 
larger effect, as shown recently by Smith and 
Keefer[21], but these effects have been neglected in 
the current model. 

The finite difference solution procedure employed 
in this work is an adaptation of the method described 
in detail by Patankar[22], A detailed description of 
the numerical procedure used in the current calcu- 
lation was given in[3]. The SIMPLE algorithm[22,23], 
which uses a conservative form of finite difference 
equation, primitive variables and staggered grids, was 
adopted to solve the coupled mass, momentum and 
energy conservation equations (1)— (4). 

3.2. Numerical results 

The validity of the present model was established 
by calculations for argon plasmas having a wide 
range of gas inlet velocities, gas static pressures and 
laser power, using two different lenses of 203 mm 
(8-in.) and 305 mm (12-in.) focal lengths. The flow 
conditions were varied over a range of static pressure 
from 1.3 to 4.0 atm, flow inlet velocity from 0.4 to 
4.5 m/s and laser power from 261 to 967 W. The wall 
and gas inlet temperatures were assumed to be con- 
stant at 300 K, and the inlet velocity was assumed to 
be uniform for these calculations. The C0 2 laser 
operated at a wavelength 10.6 pm and the collimated 
laser intensity distribution was that used by J eng and 
Keefer[3]. It is an annular beam with peak intensity 
at a radius of 25.5 mm. The dimensions of the 
203 mm (8-in.) focal length lens are given in Fig. 7. 

The semi-two-dimensional model of Glumb and 
Krier[20] assumes zero radial velocity which is equiv- 
alent to a constant axial mass flux assumption. 
However, the large temperature gradients within the 
LSP result in large density gradients and a very 
complex flow in the plasma. Figure 8 shows a typical 
calculated plasma from our two-dimensional model. 



axial distance (mm) 

Fig. 8. Isothermal contours and flux vectors for a 
calculated plasma. The contour interval is 1000 K and the 
outer contour has a value of 1000 K. 
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Fig. 9. Comparison of predicted and measured temperature 
contours for a high-velocity plasma. The contour interval is 
1000 K and the outer contour has a value of 10,500 K. 

The isotherms are shown in the top half of the figure, 
and the mass flux vectors are shown in the bottom 
half. The laser focal zone is located 50 mm from the 
inlet, and the direction and magnitude of the local 
mass flux, at the vector origins, is proportional to the 
direction and the length of the vectors. 

The LSP is positioned near the laser focal zone, 
and the peak temperature is 16,390 K. In the high 
temperature region, where the axial velocities were 
much larger than elsewhere, the axial mass flux was 
very small compared to the inlet mass flux. It was 
found that a local high pressure zone was produced 
near the leading edge of the plasma that forced the 
flow to have outward radial velocities near the center- 
line. From this result, it is clear that the constant axial 
mass flux assumption is inadequate. 

The model has been compared with a wide variety 
of experimental results given by Welle et o/.[4J in 
order to learn the capability and limitations of the 
model. Since the spectroscopic measurements were 
limited to temperatures greater than 10,000 K, Figs 9 
and 10 contain only those isotherms within that 
temperature range, and a limited portion of the 
calculation domain is plotted in order to gain a 
clearer comparison. 

The performance of the model for two different 
inlet velocities, using the 203 mm (8-in.) focal length 
lens, is demonstrated in Figs 9 and 10. The static 
pressure and incident laser power of the two cases 
differ by less than 1%, but the inlet velocity of the 
LSP in Fig. 9 is about seven times that in Fig. 10. The 
model calculates the size and position of both plas- 
mas well. As the inlet velocity increases, the radius of 
the plasma decreases, and the plasma stabilizes fur- 
ther downstream. The model predicts two local tem- 
perature maxima within the low-velocity LSP. One is 
on the laser path upstream of the focal zone, and 
another is at the focal zone. Both calculated max- 
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axial distance (mm) 


Fig. 10. Comparison of predicted and measured tem- 
perature contours for a low-velocity plasma. The contour 
interval is 1000 K and the outer contour has a value of 
10,500 K. 

imum temperatures are about 1000 K higher than the 
measured peak temperature (14,500 K). Only a single 
maximum temperature was measured, and it is lo- 
cated neither on the laser path nor at the focal zone. 

The predicted local temperature maximum near the 
focal zone, which was not observed in the experiment, 
may result from the neglect of refraction of the 
laser beam. Although it has been demonstrated 
experimentally[2] that this assumption leads to only 
a small effect on the measured global laser power 
absorption, its local effect on the predictive model 
can be more serious, since the local temperature is 
strongly influenced by the local laser intensity. When 
refraction of the laser beam within the LSP is consid- 
ered, it results in a larger focal spot size and reduced 
laser intensity. The spot size for the refracted beam 
is highly dependent on the LSP position and size. For 
a low-velocity LSP, most of the plasma is located 
upstream of the focal zone, and this makes the focal 
spot larger. When diffraction effects are combined 
with lens aberrations, the diffraction limited focal 
spot is approximately twice the size of that calculated 
from a geometric ray trace[21]. The combination of 
laser beam refraction and diffraction, and the longer 
absorption pathlength ahead of the focal zone, results 
in an actual laser intensity which is smaller near the 
laser focal zone than that which is predicted by the 
ray tracing model. Thus, the local temperature max- 
imum will not occur as predicted in the low-velocity 
LSP. This may explain why the calculations do not 
agree with the experiments in this detail. 

Since the size and temperature of these plasmas 
were well predicted, it is not surprising that the 
calculated thermal radiation power loss (316 and 
445 W for high and low velocity, respectively) agreed 
with the experimental results (326 and 392 W for high 
and low velocity, respectively). The model prediction 


for the transmitted laser power in the high-velocity 
LSP also agrees reasonably well (110 and 159 W for 
prediction and experiment, respectively). However, 
the predicted transmitted laser power (91 W) in the 
low-velocity LSP is much smaller than the mea- 
surement (246 W). This is because the predicted 
second local maximum temperature near the focal 
zone absorbs an unrealistically large amount of addi- 
tional laser power. 

In general, the model predicted well the experi- 
mentally observed changes in LSP characteristics 
resulting from changes in pressure, optical geometry 
and incident laser power for the low power argon 
plasmas, as discussed at length by Jeng et a/.[4]. Since 
the model was successful in predicting the behavior of 
low power argon plasmas, it was extended to permit 
the prediction of the higher power hydrogen plasmas 
that are of interest for applications in propulsion[5]. 
The model was modified to use the thermodynamic 
and transport properties of hydrogen and it was 
assumed that the laser beam incident on the lens had 
a Gaussian profile whose diameter could be varied. A 
parametric study was performed over the following 
range of flow and optical conditions: static pressure 
of the flow, 1-10 atm; flow inlet velocity, 5-100 m/s; 
laser power, 5-60 kW, and Gaussian laser beam 
diameter, 40-80 mm. Many examples of the model 
predictions can be found in[5], and a representative 
case is shown in Fig. 11. 

Temperature contour lines range from 1000 to 
16,000 K in 1000K increments, and the maximum 
temperature of 16,804 K. is located at x =37.4 mm, 
which is well upstream of the laser focal point 
( x = 52 mm). The magnitude of the local mass flux, 
at the vector origins, is proportional to the vector 
lengths, and the flow direction is the same as the 
vector directions. Figure 12 is a close-up view of the 
same plasma where the beam boundary is also plot- 



Fig. 11. Isothermal contours and mass flux vectors for a 
calculated hydrogen plasma. 


a aiG'.NAL r/' B 

OF POOR QUALITY B -9 




Laser thermal propulsion using laser sustained plasmas 


375 



Fig. 12. Close-up view of the isotherms and the velocity 
vectors for the hydrogen plasma of Fig. 11. 


ted. In the bottom half of this figure, the vectors 
represent the velocity vectors instead of mass flux 
vectors. Near the center axis, upstream of the max- 
imum temperature, the directions of energy transfer 
by convection and conduction are opposite, so the 
axial temperature gradient is large in this region. On 
the other hand, the temperature gradient downstream 
of the maximum temperature, where the mixing 
between inner hot and outer cold flow occurs and 
laser heat addition is small, is much smaller than that 
in the upstream region. 

The temperature ratio between the core region of 
the LSP and the inlet is as large as 50, the specific 
volume ratio can reach 150, and volumetric expan- 
sion is the dominant factor in the flow acceleration. 
In Fig. 12, the axial velocity increases from 20 m/s at 
the inlet to 255 m/s in the plasma core region, but the 



ini«t velocity (m/s) 

Fig. 13. Fractions of the incident laser power that are 
absorbed by and radiated from the hydrogen LSP as a 
function of inlet velocity. The incident laser power is 20 kW, 
the plasma pressure is 3 atm, and the Gaussian beam 
diameter is 40 mm. 


axial mass flux (see Fig. 11) is smaller in the high 
temperature region than at the inlet. Near the center- 
line upstream of the plasma, the axial mass flux 
decreases along the axial direction, so a small out- 
ward radial velocity develops (for example at 
jc*35mm and r* 3.0 mm in Fig. 12). The axial 
mass flux within the high temperature region 
(T > 6000 K) is less than one-fifth of that at the inlet, 
and this indicates that the constant axial mass flux 
assumption used by quasi-two-dimensional models is 
not adequate. 

The fractional power absorption and thermal radi- 
ation loss are shown as a function of inlet velocity in 
Fig. 13. The LSP absorbs most of the incident laser 
power and only a small fraction of laser power is 
transmitted (maximum fractional transmission is 
1.7% at the lowest inlet velocity of 5 m/s). It was 
found that the fractional power absorption increases 
as the inlet velocity increases. This is because the 
plasma has a longer tail in the high-velocity cases 
which increases the total absorption pathlength. The 
fractional power lost from the plasma through ther- 
mal radiation decreases from 72% at 5 m/s to 34% at 
lOOm/s. For optically thin radiation, the radiated 
power is proportional to the volume of the plasma, 
and since the diameter of the high velocity plasma is 
smaller, the radiation power loss decreases with in- 
creasing velocity, as shown in Fig. 13. The difference 
between power absorbed and power radiated repre- 
sents the net power converted from the laser directly 
to the flow. By increasing the inlet velocity, the net 
direct power conversion efficiency increases from 26 
to 64%. 

These results suggest that hydrogen plasmas sus- 
tained at higher laser powers with Gaussian beam 
profiles are capable of absorbing essentially all of the 
incident laser radiation, and that the thermal radi- 
ation losses can be limited to values that can be 
absorbed through the chamber walls in a regenerative 
cycle. 

4. SUMMARY 

Spatially detailed measurements have been ob- 
tained for a large number of plasmas sustained by a 
carbon dioxide laser in forced convection flows of 
argon. Analysis of these experiments have revealed 
the complex interactions between the flow, pressure 
and laser power that determine the stable operating 
regimes for these plasmas. It was found that the 
position of the plasma within the focal volume of the 
laser beam could be controlled with the incident flow 
velocity. Controlling the flow, together with the 
focusing geometry of the laser beam, permitted stable 
operation of the plasma with high fractional absorp- 
tion of the incident laser beam and relatively small 
thermal radiation losses from the plasma. A numer- 
ical model that incorporates the complex energy 
conversion processes within the plasma has been 
developed and verified using the detailed experi- 
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mental measurements. This model was then extended 
to calculate hydrogen plasmas that were sustained by 
a higher power laser having a Gaussian beam profile. 
The results of these calculations indicate that hydro- 
gen plasmas can be sustained that will absorb essen- 
tially all of the incident laser power and that lose only 
a small fraction of the absorbed power through 
thermal radiation. 

The numerical studies for hydrogen plasmas indi- 
cate that, when higher power Gaussian laser beams 
are employed, the chamber walls can be adequately 
cooled by the propellent flow in an efficient regen- 
erative cycle. No “fatal flaws” have been discovered 
by either the experimental or the theoretical in- 
vestigations, and it appears that a practical, high 
specific impulse propulsion system could be based on 
the laser thermal propulsion concept if lasers having 
adequate power are developed and if suitable window 
materials can be found. 
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ABSTRACT 

A numerical code has been successfully devel- 
oped for the investigation of thruster performance 
using a laser-sustained hydrogen plasma as the pro- 
pellant. The plasma was sustained using a 10.6/im 
CO2 laser beam which is focused at different posi- 
tions within the thruster. The physical model as- 
sumed that plasma is in thermodynamically equilib- 
rium (LTE), and geometric ray tracing was adopted 
to describe the laser beam. The steady-state, ax- 
isymmetric, Navier-Stokes equations coupled with the 
laser power absorption process have been solved nu- 
merically. A pressure based Navier-Stokes numeri- 
cal solver using body-fitted coordinates was used to 
calculate the laser-supported rocket flow which in- 
cludes both subsonic and supersonic flow regions. 
From the limited parametric study, which did not 
try to optimize the rocket performance, it was found 
that better performance was obtained when the 
laser beam was focused closer to the rocket throat. 


NOMENCLATURE 


A 

I 

J 

k 

P 

9 i» 92 , 93 
9 


coefficient in difference equation for 

variable <f> 

laser ray intensity 

Jacobian of inverse coordinate 

transformation 

conductivity 

pressure 

geometric relations between coordinate 
systems * 

emission coefficient for optically thin 
radiation 
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J2, 5, T source term in equations 6, 7, 15, 
respectively 

3 distance along laser ray 

U , V velocity component along f — and 

17— axis, respectively 

u, t; velocity component along x— and 

r— axis, respectively 

f , rj axes of curvilinear coordinate system 

<f> general dependent variable 

p density 

T effective diffusion coefficient 

/i viscosity 

a absorption coefficient 

Subscript 

E, Wy Ny Sfour adjacent nodes to P 
ty w 9 n, s four surfaces of the control volume 
centered at P 
i laser ray i 

P nodal point to be solved in difference 

equation 
rad radiation 

f , rj partial derivative with respect to f , rj 

<f> for dependent variable <j> 

Superscript 

* values based on given pressure field 

' corrected values according to corrected 

pressure field 

INTRODUCTION 


This paper describes the numerical calculations 
for a rocket using a laser-sustained plasma as the pro- 
pellant for a propulsion device. The plasma is used 
to absorb and convert the power from a remote laser 
beam into the propellant which is expanded through 
a nozzle to generate the required thrust for a space 
vehicle. This kind of propulsion system is called laser- 
thermal propulsion. Several different modes of oper- 
ating a laser-supported rocket have been proposed; for 
example, pulsed laser-supported detonation wave; 1 
cw laser-sustained plasma (LSP). 2 ” 4 In this paper, 
the discussion is limited to the numerical calculations 
of a rocket flow using <* cw laser-sustained pure hy- 
drogen plasma as the propellant, 

A schematic of a proposed laser-supported rocket 
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thruster is shown in Figure 1. Laser power is ab- 
sorbed by and sustains a high temperature plasma 
core (about 15,000-20,000 K), which mixes with the 
cold buffer flow of propellant downstream, heating the 
flow to the maximum temperature that the chamber 
wall material can tolerate. Then, the flow is accel- 
erated through a conventional nozzle to produce the 
required thrust. Since the laser power absorption by a 
plasma is through inverse bremsstrahlung, any type of 
propellant can be used, and for propulsion purposes, 
the lowest molecular weight propellant, hydrogen, is 
the best candidate. The major advantages of *bis de- 
vice, compared to chemical propulsion, are that the 
lowest molecular weight propellant can be used, and 
that the maximum temperature of the working fluid is 
not limited to the flame temperature of a combustion 
process. The combination of high temperature and a 
low molecular weight propellant results in a very high 
expected specific impulse: An estimated specific im- 
pulse for laser-sustained plasma propulsion using hy- 
drogen propellant is in the range of 1000-1500 sec., 6 
which is about three times that of the most advanced 
chemical propulsion system. 

The modeling of the LSP for propulsion began 
with the one-dimensional work of Raizer. 6 Several 
other studies 7 ’ 8 followed this analysis but used more 
sophisticated analytical and computational methods. 
Unfortunately, one-dimensional models have an inher- 
ent weakness. The LSP tends to be highly dependent 
on laser beam geometry, radial convection and radial 
diffusion transport of momentum and energy in the 
flow. Keefer and coworkers 9,10 developed a simplified 
semi-two-dimensional model which considered the ra- 
dial energy diffusion but assumed that the laser beam 
and the streamlines were parallel. Later, dumb and 
Krier 11 extended that model by using a converging 
laser beam, as well as using more realistic gas prop- 
erties. However, their model still did not consider the 
momentum equations but, instead, assumed that the 
radial velocity was zero. 

Merkle and coworkers 12 ’ 13 used a contemporary 
numerical method to attack the full two-dimensional 
problem. They succeeded in the calculation of a low 
temperature LSP (peak temperature about 4,000 K) 
in a low-mach number hydrogen flow seeded with ce- 
sium. It is feasible to extend their approach to a pure 
hydrogen plasma, but the higher peak temperature 
(about 16,000 K) and much stronger nonlinear opti- 
cal properties of the pure hydrogen plasma can result 
in numerical stability problems. 

Recently, Jeng and Keefer 14 ’ 15 obtained the 
numerical solutions for a rigorous two-dimensional 
model with realistic gas properties for the LSP in 
pure argon and pure hydrogen flows. Their results for 
the argon flows agree well with experimental data 10,17 
in the prediction of plasma position, size, shape and 


peak temperature. Jeng et al 16,18,19 have also stud- 
ied the effect of laser power, flow rate, static pressure 
and optical geometry on the LSP flows. Their work 
also indicated that there are no fatal errors in the 
laser propulsion concept, and the critical parameters 
of plasma size, position, radiation heat loss and en- 
ergy conversion efficiency, can be controlled using an 
appropriate optical arrangement and a suitable flow 
configuration. 

The objectives of the current study are to de- 
velop a new numerical code based on extended plasma 
model developed by Jeng and Keefer 14,15 for the real 
nozzle flow which includes both low and high mach 
number flows, and to conduct a preliminary paramet- 
ric study of laser rocket performance. There was no 
attempt, in this paper, to determine the actual rocket 
performance. In the following section, the theoretical 
method used in the study is described. This is fol- 
lowed by a discussion of model calculations and some 
initial parametric studies. 

THEORETICAL MODEL 

LASER BEAM. Geometric optics were used to 
describe the laser beam which was assumed to be di- 
vided into a finite number of individual rays. Each 
ray path through the optical system to the plasma 
flow was calculated using real-ray tracing. Diffrac- 
tion of the laser beam due to the finite aperture of the 
lens, and effects due to the inhomogeneous refractive 
index within the plasma were neglected in the calcu- 
lation. The local intensity of each individual ray » was 
described by Beer’s law written in the the following 
form. 


where is the distance along the laser propagation 
direction for ray t, and a is absorption coefficient at 
laser wavelength and depends on local plasma tem- 
perature and pressure. 

PLASMA FLOW. The model assumed that the 
flow is laminar, steady-state, axisymmetric and has 
variable physical properties. The pressure of the 
plasma is relatively high so local thermodynamic equi- 
librium (LTE) for the plasma can be assumed with 
little error. Therefore, the plasma can be described 
by a single temperature, and its physical properties 
are only a function of this temperature and pres- 
sure. The model also considered the thermal radia- 
tion emitted from the plasma which included the con- 
tributions from free to free, free to bound, and bound 
to bound transitions. These radiation properties are a 
strong function of wavelength, and the current model 
adopted the simplified approach, developed by Kemp 
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et al 20 and Caledonia et aL 21 to solve the hydro- 
gen plasma radiation transport equations. Radiation 
transport was divided into two parts depending on 
wavelength: for wavelengths longer than 0.095 /zm, 
the optical depth is long and an optically thin ap- 
proximation, in which all radiation escapes from the 
plasma, was used; and in other wavelength regions, 
the optical depth is short and an optically thick ap- 
proximation, which can be described as a diffusion 
process, was assumed. Following the above assump- 
tions, the equations of conservation of mass, momen- 
tum and energy for the flow can be written as : 


3(pu) 1 d(rpv) 

dx r dr 


( 2 ) 


d(puu) 1 d(rpuv) 
dx + r dr 


■ 


d ( du\ 

toSfai) 

1 d ( du\ d ( du\ 13/ dv\ 

"rTr^d;) + Tx^Tx) + -rTA^Tx) 

d [2a» (du I3(ru) \l 
dx [3 V3x r dr )y 


(3) 


where k e j / is the sum of molecular thermal conduc- 
tivity and radiation-induced thermal conductivity in 
the optically thick limit; q ra d is radiative heat loss 
in the optically thin limit; a£ 1$ represents the local 
power addition from the laser beam to the flow; and 
as mentioned before, the physical properties involved 
in equations (2-5) are a function of the local temper- 
atures and pressures. Reference 15 has a detailed de- 
scription of the thermophysical and optical properties 
sources used in this work. 

In the calculations, adiabatic, nonslip and zero 
emissivity conditions were assigned at the thruster 
wall. Total pressure, total temperature and zero ra- 
dial velocity were specified at the inlet, and upstream 
running Riemann invariants extrapolation plus the 
isentropic relations were used to calculate other prop- 
erties. At the exit, the properties were linearly ex- 
trapolated from the interior flow domain. 

METHOD OF COMPUTATION 

COORDINATE TRANSFORMATION. The set 
of transport equations (3-5) can be written in the 
cylindrical system coordinates for the dependent vari- 
able, <f> t in the following generalized form: 


dv\ 

l dx) 


d(puv) 1 d(rpvv) _ dp d / 
dx + r dr dr + dx \ 

Id/ dv\ d ( du\ Id/ dv\ 2 p 
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+r ■ R+{x, r), (6) 

where V^ e /j is the effective diffusion coefficient and 
iZ^(x, r) is the source term. The source term for the 
momentum equations (3,4) also includes the diffu- 
sion terms from the dilation effect. When new in- 
dependent variables f and rj are introduced, Eq. (6) 
changes according to the general transformation f and 
rj. A schematic illustration of the relations between 
the physical domain and the transformed domain is 
shown in Figure 2. Equation (6) is rewritten in $ and 
rf coordinates as follows: 


+rS+{£ t r,) (7) 

where 

U -ur n - vx n , 

V a* VX{ — ur^, 
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( 8 ) 

(9) 


7 1 = 3% +rjj, 

(10) 

92 = + r € r„, 

(11) 

93 = +y\, 

(12) 

J — X (r n x r) r (> 

(13) 


and S*(f, rj) is the source term in the f and rj coor- 
dinates. 

NUMERICAL ALGORITHM A nonstaggered 
grid system, as shown in Figure 2, is adopted in the 
calculations. For a typical grid node P are enclosed 
in its cell and surrounded by its neighbors N, S, E, 
and W. The integral of Eq. 7 over its control volume, 
by arbitrary taking Af = 1 and A rj = 1, becomes 


{rpU<f>) e - [rpU<f>) w + [rpV <f>) n - {rpV<f>). 

[ rT + 1 r r r* 

[ r r^ *1 r rr^ 

~ (“ 72 ^< + 73 ^,) 

+r ■ J ■ 


§ 

(14) 


Second-order center differencing is used to dis- 
cretize the above equation, except for the convective 
terms when the cell Reynolds number is greater than 
2. In that case, first order upwind differencing is used 
for the convective terms. A staggered grid system, for 
body-fitted coordinates, has been widely used, 22 and 
it has been known that an oscillatory field will be pro- 
duced using the nonstaggered grid system. The arti- 
ficial viscosity term proposed by Rhie and Chow 23 
to cure the problem was also implemented in this 
study. The cross derivatives in the diffusion terms 
of Eq. (14) are the results of the nonorthogonal coor- 
dinate system and are usually small. In order to use 
the solution procedure for the five-point unknowns, 
these terms are combined with the source term and 
treated as known quantities and the resulting rela- 
tion between <f>p and the neighboring values can be 
written as : 


A p <j)p = Ae4>e + Aw<t>w + As<I>n + As<t>s 

+rJT+{s, r,), (15) 

where the coefficients A involve the flow properties of 
convection, diffusion, area, etc., and includes 

the source term f?) and the cross derivatives in 

the diffusion terms. 

PRESSURE CORRECTION PROCEDURE. 
While the general transport equation, Eq. (15), was 


used to solved momentum and energy equations, the 
pressure field must be solved using the continuity 
equation. The procedure used in this work is an adop- 
tion of algorithm developed by Rhie, 24 and only a 
brief discussion is presented in this section. The pres- 
sure equation used in the adopted procedure is similar 
to that of SIMPLER 25 and PISO, 26,27 except that the 
density term is treated implicitly for the compressible 
flow. The basic idea is to use the following relations 
to correct the mass flux obtained from momentum 
equations 


pu = p* u* + p 9 u* + p*u' (16) 

pv = p* v* + p* v* + p* v* (17) 

and the discretized pressure equation becomes a 
transport equation for compressible flow instead of 
Poisson equation. The ratio of diffusion term to con- 
vection term is dependent on Mach number, and as 
Mach number increases, the convection term becomes 
dominant. The split operator concept similar to PISO 
was also adopted in the current study to solve the 
pressure correction equation. 

RESULTS AND DISUCSSION 

Figure 3 shows the thruster dimensions used in 
the calculations. It includes a cylindrical section up- 
stream of the convergent-divergent nozzle. The throat 
and cylindrical diameters are 6 and 40.74 mm respec- 
tively, and the exit to throat area ratio is 6. The coor- 
dinate origin used to present the calculated results is 
at the center of the throat. Although Figure 1 shows 
that the hydrogen flows into the chamber radially, 
in the calculations reported here the inlet flow direc- 
tion was assumed axial. The total pressure and total 
temperature were specified at the inlet. In the tradi- 
tional chemical rocket, the total pressure has a minor 
effect on the expected value of the specific impulse, 
but in the laser-supported rocket, the pressure has a 
substantial influence on the absorption and radiation 
processes that directly affect the available enthalpy of 
the propellant and the specific impulse of the thruster. 
Based on our previous parametric studies 15 for laser- 
sustained hydrogen plasmas in a constant area pipe 
flow, it was found that 3 atm total pressure was ad- 
equate for the thruster upstream condition, and this 
value was used in the calculations. The total tem- 
perature at the inlet was specified as 1500K in the 
calculations, based on the expected use of a regenera- 
tive cycle to recover the optically thin radiation from 
plasma. 

The focusing lens used in the calculations has a 
nominal focal length 203 mm, and References 14 and 
15 has detailed discussion of this lens geometry. The 
actual dimensions of this lens were adopted in the ray 
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tracing for the laser beam, so the spherical aberration 
of the lens was also included in the calculations. The 
laser was assumed to have 10.6 /zm wavelength and a 
collimated, 40 mm diameter Gaussian beam incident 
on the lens. The only preliminary parametric study 
of the rocket performance reported in this paper was 
to focus the laser beam at different axial locations 
using the same lens arrangement. Table 1 lists the 
operating conditions, including one case without laser 
power addition. 

In the calculations, nonuniformly spaced grids in 
the physical space were transformed into uniformly 
spaced grids in both radial and axial directions in the 
transformed space, as illustrated in Figure 4. From 
previous studies 14 »15*1S*19 Q f LSP within a constant 
diameter pipe, the position of the high temperature 
plasma zone was found to be highly dependent on the 
flow configurations and the optical arrangements, and 
it was not easy to prelocate the region in the physical 
plane where a dense grid distribution was needed to 
accommodate the high temperature gradients. There- 
fore, a preliminary calculation was first performed for 
each case to locate the high temperature zone, then 
the grid distribution was refined and used for the final 
calculation. 

The newly developed computer code was first 
used to calculate the thruster flow without consid- 
ering laser heat addition, and the results were used 
to check the behavior of the computer program and 
to serve as benchmark to compare with the other 
calculated cases which included the laser heat ad- 
dition. The Mach number contours for case no. 1 
are shown in Figure 5. Because the nozzle expan- 
sion angle is only 6°, the contours outside of the wall 
where diffusion transport is not important, are almost 
straight, and are similar to one-dimensional nozzle 
calculations. Since the nonslip condition was applied 
in the calculations, the Mach number near the wall is 
smaller than that near the centerline at the same axial 
location. This demonstrates that the newly developed 
computer program can also handle the boundary layer 
in supersonic flow. In this test case, there were only a 
few grids within the boundary, and we made no spe- 
cial attempt to resolve the boundary layer to good 
numerical accuracy. If needed, more grids can be in- 
cluded in that region and a more accurate solution 
can be obtained. 

Figure 6 shows isotherms for case no. 2, and the 
outer laser beam boundary was also plotted in the 
figure as a dashed line. The plasma core is located 
upstream of the laser focus, and the distance between 
the laser focus and the high temperature plasma core, 
which depends on flow configurations and optical ar- 
rangements as discussed in the early section of this 
paper, is 10 mm. Downstream of the laser focus, the 
mixing process quickly reduces the centerline temper- 


ature. In the throat region, only a small portion of 
the hydrogen flow near the axis is heated, and a major 
portion of the flow is still unheated. For an optimized 
rocket, one of the design criteria is to have a uniformly 
distributed temperature at the throat in order to re- 
duce the two-dimensional losses. In the calculated 
case, this condition was not satisfied. Near the wall 
of the nozzle exit, there is a local hot spot which is due 
to the viscous dissipation effect within the boundary 
layer. 

Laser-sustained plasmas, operated at 30 k W laser 
power, absorb almost all of the incoming laser beam 
power. In the calculated case no. 2, only 397 W 
of laser power was transmitted through the plasma 
which either reached the thruster wall or was trans- 
mitted through the throat. The wall material could 
easily withstand this small level of transmitted power 
density. The optically thin radiation from the plasma 
is also an important factor for a laser-supported 
rocket design. It represents the portion of laser power 
that is not directly converted into the flow power. 
However, this radiation power, like a part of the trans- 
mitted laser power, will be absorbed by the wall, and 
by using a regenerative cycle, part of this energy can 
be recovered and used to heat the hydrogen from its 
storage temperature up to the propellant inlet tem- 
perature of the thruster. For case no. 2, the optically 
thin radiation was 17.9 kW which was only 30 % of 
the energy needed to heat hydrogen to the assigned 
inlet temperature, 1500 K. 

Pressure and Mach number contours for case no. 

2 are shown in Figure. 7 and Figure. 8, respectively. 
Pressure contours downstream of the laser focus in 
Figure 7 are similar to the calculated results for case 
no. 1. Those have only a very small radial variation 
and imply that the nonuniformly distributed temper- 
ature field has little effect on the pressure field. How- 
ever, a high pressure zone, even higher than the inlet 
static pressure of the thruster (3xl0 5 N/m 2 ) } occurs 
ahead the high temperature plasma zone. The Mach 
number at the inlet is about 0.019 and the maximum 
calculated Mach number is 2.605 at the nozzle exit. 
The Mach number contours in the supersonic region 
are very complicated, compared to that of case. no. 

1 in Figure 5, although the velocity monotonically in- 
creased from the wall to the centerline. The Mach 
number at the centerline is always smaller than away 
from the centerline, since the temperature, as well as 
sonic speed, is higher at the centerline. The unheated 
region still has very straight Mach number contours 
similar to those shown in Figure 5 for case no. 1. The 
mixing region between hot center flow and the colder 
unheated flow shows very complicated contours. Mix- 
ing in that region, which implies variations of sonic 
speed due to different temperature, molecular weight, 
and specific heat, are the major factors. The Mach 
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number near the wall is smallest due to the bound- 
ary effect, and the Mach number on the centerline is 
always smaller than that away from the centerline in 
the subsonic portion of the flow. 

In the previous studies 15,16 , it was found that 
laser heat addition to the flow will generate a very 
complicated velocity field, even in a constant area 
pipe. The same behavior was also observed in the 
current results. Figure 9 shows the radial velocity 
contours for case no. 2. Near the high temperature 
plasma core region, both positive and negative radial 
velocity occurs. In the supersonic region, a local min- 
imum exists In the mixing zone between heated and 
unheated flow for case no. 2, and the radial veloc- 
ity increases from the centerline up to the edge of 
the boundary for the unheated case no. 1. Nonuni- 
formly distributed temperature at the throat is the 
main cause for the above results. 

Figure 10(a) and 10(b) show isotherms for cases 
no.4 and no. 5, respectively. Their laser focuses are 
closer to the throat than case no. 2, and the focus is 
within the nozzle contracting section for case no. 5. 
As the laser focus moves closer to the throat, the size 
of the plasma becomes smaller. Since the optically 
thin radiation is proportional to the plasma size, the 
radiation loss for case no. 5, 9.9 kW, is much smaller 
than 15.9 kW and 18.6 kW for cases no.4 and no. 
2, respectively. Because the transmitted laser power 
is small and similar for all cases, the energy directly 
converted to the flow is highest for case no. 5. The 
laser beam also heats a larger portion of the flow at 
the throat in case no. 5 than other cases, and this 
produces better thruster performance. 

Table 2 lists the calculated properties for the 
five test cases. The optically thin radiation heat loss 
ranges from 61 % to 33 % of the laser power, and the 
transmitted laser power is relatively small by com- 
parison. The peak axial velocity, which is always lo- 
cated on the centerline at the exit, increases from 4763 
m/sec for the unheated thruster to 18957 m/sec for 
case no. 5. The maximum temperature of the plasma 
shows little dependence on the variation of laser focal 
position among the test cases. The peak Mach num- 
ber at the exit plane is located within the unheated 
flow for all test cases, and it increases as the laser 
focal point moves downstream. The mass flow rate 
decreases as the laser focal point moves toward the 
throat. 

SUMMARY AND CONCLUSIONS 

The flow within a laser-sustained plasma rocket 
thruster has been successfully calculated. The physi- 
cal model was based on steady-state two-dimensional 
Navier-Stokes equations for the flow and geometric 
real-ray tracing for the laser beam. Realistic thermo- 


physic and optical properties based on locally ther- 
modynamical equilibrium ^sumptions were incorpo- 
rated into the calculations, and radiation heat trans- 
fer from the plasma was simplified using optically thin 
and optically thick assumptions. The governing equa- 
tions were first transformed into body-fitted coordi- 
nates and were solved numerically. A pressure-based 
Navier-Stokes solver with a split operator concept was 
used for the numerical procedure. 

Several laser propulsion cases were calculated us- 
ing the newly developed computer code. The noz- 
zle geometry, flow configuration and optical, arrange- 
ments were not optimized for rocket performance. 
However, some useful information was obtained from 
these preliminary test cases. The performance of the 
rocket is highly dependent on the laser focal position. 
The maximum velocity is larger when laser focal point 
moves closer to the throat. Generally speaking, the 
laser beam only heated a small portion of the hydro- 
gen propellant, and the best thruster performance will 
not be achieved due to the two-dimensional losses. In 
the future, an extensive parametric study will be re- 
quired to optimize the rocket performance, and one 
of the major goals of the study is to design the rocket 
to produce a more uniformly distributed temperature 
at the throat. 
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Figure 10. Temperature contours: (a) case 
number 4, (b) case number 5. 


Table 2: Calculated properties for five test cases 


Case no. 


Transmitted power (W) 
Max. axial vel. (m/sec) 
Max. temperature (K) 
Max. mach number 
Mass flow rate (v/sec) 
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ABSTRACT 

An extensive numerical experiment has been con- 
ducted to evaluate rocket thruster performance us- 
ing a laser-sustained hydrogen plasma as the propel- 
lant. The plasma was sustained using a 30 kW CO 2 
laser beam operated at 10.6 fim focused inside the 
thruster. The steady-state Navier-Stokes equations 
coupled with the laser power absorption process have 
been solved numerically. A pressure based Navier- 
Stokes solver using body-fitted coordinate was used 
to calculate the laser-supported rocket flow which in- 
cluded both recirculating and transonic flow regions. 
The local thermodynamic equilibrium (LTE) assump- 
tion was used for the plasma thermophysical and opti- 
cal properties. Geometric ray tracing was adopted to 
describe the laser beam. Several different throat size 
thrusters operated at 150 and 300 kPa chamber stag- 
nation pressure were studied. It was found that the 
thruster performance (vacuum specific impulse) was 
highly dependent on the operating conditions, and a 
properly designed laser supported thruster can attain 
a specific impulse around 1500 secs. The heat loading 
on the thruster wall was also estimated and was in the 
range of that for a conventional chemical rocket. 

NOMENCLATURE 

I local laser ray intensity 

k conductivity 

p pressure 

q optically thin radiation 
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s distance along laser ray 

u, t; velocity component along x— and 

r— axis, respectively 
p density 

a absorption coefficient 

Subscript 

% for laser ray t 

rad radiation 

INTRODUCTION 

A propulsion system supported by a laser- 
sustained plasma (LSP) has been the subject of many 
studies. A schematic of a continuously laser sup- 
ported thruster is shown in Figure 1. The high power 
laser beam from a remote site is focused into the 
thruster. The plasma is ignited and continuously sus- 
tained near the laser focus. The plasma is used to ab- 
sorb the power from the laser beam and convert it into 
propellant enthalpy. Since the focusing lens is posi- 
tioned at the end of the thruster, the propellant flows 
into chamber radially through annular slots. The con- 
cept of two separated propellant inlets are used for the 
rocket shown in the Figure 1. The propellant flowing 
through the slot near the lens is relatively cold and 
serves to cool the lens. The major propellant flow 
enters the thruster through the other slot and has a 
relatively high temperature. This temperature is de- 
pendent on the regenerative heat from the thruster 
wall to heat the propellant which was originally at 
storage temperature. Within the thruster, the flow 
convects high enthalpy plasma downstream and mixes 
with other propellant to obtain a more uniform radial 
temperature distribution. Then the propellant ther- 
mal energy is converted to kinetic energy and pro- 
duces the required thrust for the rocket through a 
converging and diverging nozzle. The major advan- 
tage of the laser thruster, compared to the chemi- 
cal thruster, is higher specific impulse (about 1000- 
2000 secs) since the lowest molecular weight propel- 
lant (hydrogen) and a higher propellant temperature 
(not limited to the flame temperature of a combustion 
process) can be used. 
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The modeling of the LSP for propulsion be- 
gan with the one-dimensional work of Raizer. 1 Jeng 
and Keefer 2 have extensively reviewed the early an- 
alytical work using one-dimensional and semi-two- 
dimensional models. A full two-dimensional model 
on the LSP was first numerically solved by Merkle 
and coworkers. 3,4 The objective of their work was to 
demonstrate the capability of the developed numeri- 
cal method, and part of the governing process of LSP 
model was simplified. For example, the plasma radia- 
tion process was neglected in their model. They suc- 
ceeded in the calculation of a low temperature LSP 
(peak temperature about 4,000 K) in a hydrogen flow 
seeded with cesium. They adopted a rocket geometry, 
but the calculations were limited to the subsonic flow 
region. It is feasible to extend their approach to a 
pure hydrogen plasma, but the higher peak temper- 
ature (about 16,000 K) and much stronger nonlinear 
optical properties of the pure hydrogen plasma can 
result in numerical stability problems. 

Jeng and Keefer 2 obtained the numerical solu- 
tions for a rigorous two-dimensional model with real- 
istic gas properties for the LSP in pure argon subsonic 
pipe flow. Their results for the argon flows agree well 
with experimental data 5,6 in the prediction of plasma 
position, size, shape and peak temperature. Jeng and 
Keefer 7,8 also extended the model for the LSP in pure 
hydrogen subsonic flows and studied the effect of laser 
power, flow rate, static pressure and optical geometry 
on the LSP structures. A study using different laser 
beam profiles and wavelengths was also investigated 
recently 0 . Their work indicated that there are no fatal 
errors in the laser propulsion concept, and the critical 
parameters, for example plasma size, position, radia- 
tion heat loss and energy conversion efficiency, can be 
controlled using an appropriate optical arrangement 
and a suitable flow configuration. 

Recently, Jeng and Keefer 10 developed a numeri- 
cal method to solve the laser supported rocket. They 
have successfully calculated a hydrogen LSP in the 
transonic flow region. They considered an axial in- 
flowing propellant, but only a limited number of rock- 
ets were investigated. Although they did not optimize 
to the rocket design, the predicted rocket performance 
is rather promising. 

The objectives of the current study are to ex- 
tend the developed numerical code for calculating re- 
alistic rocket thrusters, and to perform extensive nu- 
merical experiments for an optimized rocket design 
powered by a 30 kW C0 2 laser. The calculated flow 
structure, rocket performance (vacuum specific im- 
pulse and thrust) and the heat loading on the thruster 
wall will be discussed. In the following section, the 
theoretical method used in the study is described. 


THEORETICAL MODEL AND 
METHOD OF COMPUTATION 


An extensive discussion of the theoretical model 
and numerical approach can be found in References 
7-10. Only a brief description of the model and nu- 
merical method is presented in this section. 


LASER BEAM. Geometric optics were used to 
describe the laser beam which was assumed to consist 
of a finite number of individual laser rays. Each laser 
ray path through the optical system to the plasma 
flow was calculated using real-ray tracing. Diffraction 
of the laser beam due to the finite aperture of the lens, 
and effects due to the inhomogeneous refractive index 
within the plasma were neglected in the calculation. 
Beer’s law was used to calculate the local intensity for 
each individual ray: 


dli 


dsi 


(i) 


PLASMA PROPERTIES. Since the calculated 
rocket has a relatively high pressure, local thermo- 
dynamic equilibrium (LTE) for the plasma can be as- 
sumed with little error. Therefore, the plasma can 
be described by a single temperature, and its physi- 
cal properties are only a function of this temperature 
and pressure. Reference 7 has a detailed description 
of the thermophysical and optical properties used in 
this work. 


FLUID FLOW. It was assumed that the flow is 
steady-state, laminar and axisymmetric. The conser- 
vation of mass and Navier-Stokes equations for com- 
pressible and variable properties flow were used in the 
analysis. These equations are similar to those used in 
typical fluid mechanics analysis and are discussed in 
the previous work. 10 

One of the major interactions between a high 
power laser beam and plasma is the energy absorp- 
tion process. The energy equation, EQ, [2], used in 
the calculations is listed below, and the ^2al{ term 
represented the local heat addition from the laser to 
the fluid. 

d(puh) 1 d(rpvh) _ d fk e /fdh\ 
dx *** r dr dx \ c p dx ) 


1 d ( rk e ff dh\ _ djpuj^xP + §« 2 
r dr \ c p dr) dx 


)) 
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RESULTS AND DISUCSSION 


+ XT $ rad 

This equation also included the thermal radiation 
emitted from the plasma and included the contribu- 
tions from free to free, free to bound, and bound to 
bound transitions. These radiation properties are a 
strong function of wavelength, and the current model 
adopted the simplified approach, developed by Kemp 
et al 11 and Caledonia et al. 12 to solve the hydro- 
gen plasma radiation transport equations. Radiation 
transport was divided into two parts depending on 
wavelength: for wavelengths longer than 0.095 /xm, 
the optical depth is long, and an optically thin ap- 
proximation in which all radiation escapes from the 
plasma, was used; and in the short wavelength re- 
gion, the optical depth is short and an optically thick 
approximation, which can be described as a diffusion 
process, was assumed. The q ra d term in Equation [2] 
is radiative heat loss in the optically thin limit, and 
the k e /f is the sum of molecular thermal conductiv- 
ity and radiation-induced thermal conductivity in the 
optically thick limit. The viscous dissipation is rep- 
resented by DISS in Equation [2]. 

NUMERICAL METHOD. The set of transport 
equations for fluid flow written in the conservative 
form and in the cylindrical system coordinates was 
first transformed to the generalized curvilinear co- 
ordinate system. 10 A nonstaggered grid system is 
adopted in the calculations for the set of the transport 
equations. For the transport equations of momentum 
and energy, second-order central differencing is used 
to discretize the equations except for the convective 
terms when the cell Reynolds number is greater than 
2. In that case, first order upwind differencing is used. 
The staggered grid system in body-fitted coordinates 
has been widely used in the literature, 13 and it has 
been known that an oscillatory field will be produced 
using the nonstaggered grid system. The artificial 
viscosity term, proposed by Rhie and Chow 14 to cure 
the problem, was also implemented in this study. 

While the above numerical algorithm was used to 
solve the momentum and energy equations, the pres- 
sure field must be solved using the continuity equa- 
tion. The procedure used in this work is an adoption 
of the algorithm developed by Rhie. 15 This algorithm 
is similar to that of SIMPLE 16 and PISO, 17,18 except 
the density term is treated implicitly for compressible 
flow. The split operator concept similar to PISO was 
also adopted in the current study to solve the pressure 
correction equation. 


Figure 2 shows the thruster dimensions used in 
the calculations. As discussed in the previous section 
of this paper, the propellant temperature and flowrate 
of each inlet slot of the two-slot thruster shown in the 
Figure 1, are different and will depend on the require- 
ments for lens cooling and the regenerative heating 
efficiency. This information is not well known yet, so 
the one-slot thruster shown in Figure 2 was chosen 
for the computation. The cylindrical section of the 
thruster has a 20 mm radius, which :? the same as the 
lens and incoming laser beam radius, and 25° and 15° 
were used for the thruster converging and diverging 
angles, respectively. The other parameters used for 
the calculated rockets are tabulated in Table 1. The 
lens had a 54.5 mm nominal focal length and focused 
the incoming gaussian distributed 30 kW laser beam 
to 0.5 mm upstream of the throat. 

Nonuniformly spaced grids, 70 by 44 grids, were 
used in the calculation, and the same procedure used 
in the previous work 10 to locate a dense grid distribu- 
tion around the high temperature gradient zone was 
also implemented. The total pressure, total tempera- 
ture and flow vector were specified at the interface be- 
tween the slot and thruster chamber as the boundary 
conditions for the computations. Either 300 kPa or 
150 kPa inlet total pressure was used for the paramet- 
ric study, and zero inlet axial velocity was used. The 
total propellant temperature at the inlet is not easy to 
determine. It is influenced by the thruster wall ma- 
terial (emissivity), regenerative system performance 
and thruster internal heat transfer (convection and 
radiation). The total temperature of the inlet propel- 
lant was specified as 1500K in the calculations. This 
was based on the maximum expected use of a regen- 
erative cycle to recover the optically thin radiation 
from plasma and the practical temperature limitation 
of the thruster wall material. The thruster wall tem- 
perature was also assigned as uniformly distributed 
1500 K in the calculations. 

The calculated isotherms and normalized stream- 
lines for rocket No. 1 are shown in Figure 3 and Fig- 
ure 4, respectively. The plasma is positioned near 
the throat with a peak temperature of 15,800 K up- 
stream of the laser focus. Most of the propellant 
within the thruster is not directly heated by the laser 
beam. However, nearly all of the propellant mixes 
with the hot plasma and raises the propellant bulk 
temperature up to approximately 10,000 K before it 
enters the diverging section of the thruster. Eighty 
four percent of laser power (25 kW) is absorbed by 
the propellant, and most of the remaining laser power 
is transmitted through the throat, except for a very 
small portion that reaches the thruster wall. 
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Since the propellant enters the chamber radially 
through a slot away from the lens, part of the pro- 
pellant is expected to recirculate within the thruster. 
As shown in Figure 4, two recirculation zones occur 
in this thruster. A large one exists near the focus- 
ing lens, and the other one is near the thruster wall 
downstream of the inlet slot. The streamlines near 
the centerline upstream of the plasma indicate the 
flow has an outward direction. This is because the 
- downstream hot plasma propellant (low density) gen- 
erates a locally unfavorable axial pressure gradient in 
the flow. This phenomena was also observed in the 
previous work 2,7,8 for the LSP within a constant area 
pipe. 

An enlarged view of the isotherms and constant 
axial velocity contours around the throat are shown 
in Figure 5 and Figure 6, respectively. The leading 
edge of the LSP is 11 mm upstream of the throat. 
The axial temperature gradient at the leading edge 
is very high since the axial heat transfer of convec- 
tion and conduction is in opposite directions. The 
detailed discussion of how LSP structures are influ- 
enced by forced convection and optical geometry can 
be found in the previous work. 2 ’ 7 * 8 At the throat the 
propellant is very hot, approximately 13,500 K on the 
centerline. The temperature decreases radially from 
the centerline and reachs a value near 1250 K, which 
corresponds to the sonic temperature of the adiabatic 
(without heat addition from the laser beam) propel- 
lant. Then the temperature increases due to the wall 
heat transfer and viscous dissipation effects. 

The traditional chemical rockets are preferred to 
have a uniform temperature distribution at the throat 
to minimize the two-dimensional loss of thruster per- 
formance (Isp). For rockets supported by a laser 
beam, the spatial distribution of heat addition is 
highly nonuniform, and a nonuniform temperature 
distribution is likely to exit at the throat as shown 
in Figure 5. Actually, it is required to have some 
kind of nonuniform temperature distribution around 
the throat for a practical high performance laser sup- 
ported thruster (around 10,000 K propellant temper- 
ature at the throat), or the thruster wall will be dam- 
aged. The heat loading on the thruster wall, which 
will be discussed later in this paper, is primary due to 
irradiation from the plasma. A stream of cold propel- 
lant near the wall as shown in Figure 5, which results 
in some loss of rocket performance, is required to cool 
.. the wall. 

Upstream of the leading edge of the plasma a lo- 
cal minimum in the axial velocity occurs on the cen- 
terline where the propellant is unheated by the laser 
beam. The radial distribution of the axial velocity 
within the center plasma core region is relatively uni- 
form with the peak axial velocity on the centerline. 
Large radial variations of the axial velocity is found 


near the thruster wall where the large radial temper- 
ature gradient is also located. The exit temperature 
and axial velocity on the centerline are 10,700 K and 
25,300 m/s, respectively. The predicted vacuum spe- 
cific impulse for this rocket, with exit area expansion 
ratio of 4, is 960 secs. With further expansion of this 
rocket, the estimated specific impulse of this rocket 
can be approximately 1500 secs. 

The isotherms for two other throat radius 
thrusters operated at 150 kPa pressure, rockets Nos. 
2 and 3, are shown in Figure 7. The leading edge 
and size of these plasmas are different. The mass 
flowrate, which is directly related to the throat size, 
is the major controlling factor for the plasma size and 
position. The mass flow rate of a constant total en- 
thalpy fluid through a choked nozzle is proportional 
to the throat area. For the laser supported thruster, 
the propellant temperature and its distribution at the 
throat are other controlling factors for the mass flow 
rate. For a smaller throat thruster, a larger portion 
of the propellant at the throat has a high temper- 
ature and this will further reduce its critical mass 
flowrate. The calculated mass flow rate for cases 1-3, 
tabulated in TABLE 2, result from this effect. For 
a smaller throat thruster, the reduced axial velocity 
within the subsonic converging section of the thruster 
allows the plasma to move further upstream along the 
laser beam and to obtain a larger size, as shown in 
Figure 7. A detailed discussion of the interaction be- 
tween the plasma structure and a forced convection 
flow can be found in References 2,7 and 8. 

In the traditional chemical rocket, the total pres- 
sure has a minor effect on the expected value of spe- 
cific impulse. However, in the laser-supported rocket 
the pressure has a substantial influence on the ab- 
sorption and radiation processes which directly af- 
fect the available enthalpy of the propellant and the 
heat loading on the thruster wall. Figure 8 shows 
the isotherms for a 300 kPa inlet pressure thruster, 
rocket No 4. The absorption coefficient and optically 
thin thermal radiation for a hydrogen plasma are ap- 
proximately proportional to the 1.5 power of pressure. 
With a higher absorption coefficient in this plasma, 
compared to those of rocket Nos. 1-3, only a very 
small amount of the laser power directly escapes from 
the thruster (see Table 2); however, the optically thin 
radiation from the plasma also increases. This re- 
sults in a low conversion efficiency from laser power 
to available thermal power of the propellant. The 
higher absorption coefficient also results in a plasma 
which is positioned further upstream against the laser 
beam. Two local temperature maxima were predicted 
for this thruster. This results from two competitive 
processes which determine the local laser intensity: 
the attenuation by /p /vfill /eject /nd plasma ab- 
sorption and the increase due to the converging laser 
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beam. References 7-8 contain a detailed discussion of 
this subject. 

The vacuum specific impulse variations with area 
expansion ratio for the four calculated thrusters are 
shown in Figure 9. The most advanced chemical 
rocket has an Isp of about 500 secs, and the calculated 
Isp for the four laser supported rockets are between 
700 and 1020 secs. At the exit of these rockets the 
ratio between thermal energy and the kinetic energy 
of the propellant is still relatively high, and further 
expansion of these rockets will greatly enhance the 
rocket performance. Although the propellant within 
the expansion section is expected to be far away from 
chemical equilibrium, the reasonable Isp is expected 
to be around 1500 secs for higher expansion ratios 
using rockets Nos. 1 and 3. The highest Isp among 
the test cases is 1022 secs from rocket No. 3. As dis- 
cussed before, this thruster has the propellant almost 
fully heated at the throat. Thruster No. 1 has a lit- 
tle larger throat than rocket No. 3. and has a very 
small cold propellant stream near the throat. The 
calculated Isp of this rocket is only 20 secs less than 
thruster No. 3. With further increase of throat size, 
the size of cold stream propellant increases, and the 
calculated Isp decreases rapidly. The Isp for thruster 
No. 2 is an example. Rocket No. 4 operated at 
300 kPa and the plasma occurs further away from 
the throat than the lower pressure thrusters. This 
results in a larger cold propellant stream around the 
throat and in a much lower Isp than the lower pres- 
sure thruster. 

The heat loading on the inner wall of a good laser 
thruster should come from the plasma radiation, and 
a cold stream of propellant near the inner surface will 
reduce the heat flux through the wall. This heat flux 
is used to heat the propellant from its storage en- 
thalpy up to the inlet enthalpy. The irradiation along 
the wall of the calculated rockets are shown in Fig- 
ure 10. The irradiation first increases from the end of 
the thruster to where the hot plasma is located, and 
then decreases monotonically along the axial distance 
for all test rockets. The maximum irradiation on the 
wall for all the rockets is almost the same, although 
the 300 kPa rocket has much higher total irradiation 
flux on the wall (see Table 2). In order to estimate 
the radiation heat loading on the wall, the emissiv- 
ity and absorptivity of the wall surface are needed. 
If absorptivity is unity, the irradiation will be equal 
to the wall heat loading. Generally, the absorptivity 
is less then one, and the radiation heat loading will 
be smaller and spatially smoothed. The calculated 
maximum irradiation represents an upper limit of the 
radiation heat loading, and all maxima are much less 
then the heat loading for chemical rockets. 

Heat loading on the wall is also related to the 
thruster inlet propellant temperature which is regen- 


eratively heated from the storage temperature. The 
propellant inlet temperature will influence the rocket 
performance as long as the cold stream of propellant 
exists near the wall at the nozzle. Rockets Nos. 1 and 
3 have a very limited amount of unheated propellant 
at the throat, and the Isp for these rockets should 
not be strongly influenced by the inlet temperature. 
On the other hand, higher or lower inlet tempera- 
ture will either increase or decrease the expected Isp 
significantly for rocket Nos. 2 and 4. An inlet tem- 
perature of 1500 K was used for the test rockets, and 
the expected inlet temperatures vs. three recovery ef- 
ficiency of irradiation are tabulated in Table 3. Low 
pressure rockets Nos. 1 and 2 need the recovery ef- 
ficiency of 30 and 52 percent, respectively, for the 
adopted inlet temperature, 1500 K. If the irradiation 
is fully converted to the propellant enthalpy through 
regenerative heating for these two thrusters, the inlet 
propellant temperatures will be 2810 K and 2390 K, 
respectively. For rockets Nos. 3 and 4, the inlet tem- 
peratures will be around 3400 K, which is not realistic 
in a practical device, and probably a radiator would 
be needed to dissipate the wall heat flux and lower 
the inlet propellant temperatures. 

SUMMARY AND CONCLUSIONS 

The flow within a laser-sustained plasma rocket 
thruster has been successfully calculated based on 
a previously verified physical model 2,7 ’ 8 ’ 10 . The 
adopted numerical algorithm has demonstrated the 
capability to handle complicated flows within a 
thruster. The performance of several realistic 
thrusters powered by a 30 kW laser has been stud- 
ied. The results show that the laser supported rocket 
can have much higher vacuum specific impulse than 
the traditional chemical rocket. The Isp was in the 
range of 700 to 1020 secs for four calculated thrusters 
having an area expansion ratio of four. The irradi- 
ation from the plasma to the thruster wall was also 
calculated, and the estimated maximum heat loading 
on the thruster is within the range of the chemical 
rockets. Among four calculated rockets, rocket No. 1 
was the best practical rocket. This rocket has nearly 
the highest vacuum specific impulse among the test 
rockets, and both the estimated thruster heat loading 
and the inlet propellant temperature are reasonable. 
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Cold Hydrogen 


Hot Hydrogen 


THROAT RAOIUS- 1.2 mm 
PRESSURE- 150 k Pa 



Figure 2* Thruster geometry adopted for the com- 
putation 
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Figure 6. Axial velocity contours for rocket No. 1 
(starting from 50 m/s to 24,050 m/s with 2,000 m/s 
increments) 
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Figure 3. Isotherms for rocket No. 1 (starting from 
2,000 K to 15,000 K with 1,000 K increments) 
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Figure 4. Normalized streamlines for rocket No. 1 
(starting from -0.09 to 1.01 with 0.1 increments) 
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Figure 7. Isotherms for rocket (a) rocket No. 2, (b) 
rocket No. 3 (starting from 2,000 K to 15,000 K with 
1,000 K increments) 
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Figure 8* Isotherms for rocket No. 4 (starting from 
2,000 K to 15,000 K with 1,000 K increments) 



AREA RATIO 

Figure 9. Vacuum specific impulse vs. area expan- 
sion ratio 
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Figure 10. Irradiation along the thruster wall 


TABLE!. 

Test Conditions 


ROCKET 

1 

2 

3 

4 

Throat Radius 
(mm) 

1.2 

1.3 

1.1 

0.8 

Inlet Pressure 
(kPa) 

150 

150 

150 

300 


TABLE II. 

Performance of the Tested Rockets 


ROCKET 

1 

2 

3 

4 

Transmitted Laser 
Power (kW) 

4.9 

5.7 

4.2 

0.8 

Total Irradiation 
on Thruster (kW) 

7.1 

5.1 

11.0 

19.6 

Mass Flow Rate 
(mg / sec) 

90 

118 

73 

121 

Vacuum Specific 
Impulse* (sec) 

960 

859 

980 

686 

Thrust* 

(N) 

0.85 

0.99 

0.70 

0.80 


• at expansion area ratio of 4 


TABLE III. 

Expected Hydrogen Temperature vs. 
the Recovery Efficiency of Irradiation on Thruster 


\ROCKET 

EFPCIENcT'"^ 

r 

2* 

3* 

4** 

100% 

2810 K 

2390 K 

3310 K 

3430 K 

50% 

2290 K 

1420 K 

2830 K 

2950 K 

25% 

1295 K 

700 K 

2220 K 

2350 K 


*at150kPa 

**at300kPa 
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Laser-Sustained Plasmas in Forced Convective Argon Flow, 
Part II: Comparison of Numerical Model with Experiment 


San-Mou Jeng,* Dennis R. Keefer,! Richard Welle,! and Carroll E. Peters§ 
University of Tennessee Space Institute , Tullahoma , Tennessee 


A two-dimensional laser-sustained plasma model, which is based on the laminar, Navier-Stokes equations for 
the flow and geometric ray tracing for the laser beam, has been evaluated and compared with existing experimen- 
tal results for a wide range of forced convective argon flows. The influence of gas inlet velocity, gas pressure, 
laser power, and focusing geometry on the structure of the plasma was examined. The model agreed well with 
the existing experimental data in both global structure and detailed temperature distribution, particularly for 
static pressures greater than 2 atm. It was found that the diffusion approximation for the optically thick portion 
of the thermal radiation was not adequate for low-pressure (less than 2 atm) plasmas and that the radiation- 
induced thermal conductivity had to be adjusted in order to obtain agreement between the model calculations 
and experimental results. The present model calculations were also compared with a recently published semi- 
two-dimensional model and the results indicate that the existing one-dimensional and semi-two-dimensional 
models do not provide adequate solutions for the laser-sustained plasma. 


Nomenclature 

c p = specific heat at constant pressure, J/kg-K 
h - specific enthalpy, J/kg 

/ = laser intensity, W/m 2 

k = intrinsic thermal conductivity, W/m-K 

k iff = effective thermal conductivity, W/m-K 

k fMi = radiation-induced thermal conductivity, W/m-K 
q fMl = radiation heat loss, J/m J -s 
r = radius, m 

s = distance along the laser ray, m 

u « axial velocity, m/s 

v * radial velocity, m/s 

x = axial distance, m 

y = radial distance, m 

a = absorption coefficient at 10.6 pm wavelength. I/m 
p = viscosity, kg/m-s 

p * density, kg/m' 


Introduction 

T HE main objective of the present investigation was to 
evaluate a proposed new model* for predicting both the 
global and detailed structure of laser-sustained plasmas 
(LSP) under various flow conditions and optical geometries. 
Although the low-molecular-weight hydrogen is the best pro- 
pellant for laser thermal propulsion, there is very little 
hydrogen LSP data available at present. Therefore, this in- 
vestigation was limited to axisymmetric LSP in flowing 
argon, since we have made detailed measurements of their 
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structure in our laboratory. 2 - 3 A separate theoretical study, 
using a similar model, for the hydrogen LSP can be found in 
Ref. 4. 

It has been pointed out 5 - 6 that a basic understanding of 
the interaction between a high-power laser beam and a 
plasma is needed and that the following questions must be 
answered before designing a laser thermal propulsion 
thruster. How can the position and size of the LSP, which 
are clearly related to the size and shape of the absorption 
chamber, be controlled? How can the power conversion effi- 
ciency from laser energy to propellant thermal energy be 
maximized? How can damage to the absorption chamber 
from the high-power laser beam and the high-temperature 
plasma be avoided? In order to answer these questions, ex- 
periments with laboratory-sized LSP under various flow con- 
ditions and optical geometries have been performed. 2 - 3 
However, a theoretical model, which can accurately predict 
the properties of the laboratory-sized LSP, is needed in order 
to fill the gap between the laboratory-sized LSP and the 
larger LSP that would be required for a practical propulsion 
device. 

Since the LSP has highly nonlinear thermodynamic, 
transport, and optical properties and an extremely large ab- 
sorption of laser energy by the flow, a solution of this prob- 
lem is not easily achieved. The modeling of the LSP began 
with the one-dimensional work of Raizer. 7 Several other 
studies* -10 followed this analysis, but used more sophisti- 
cated analytical and computational methods. Unfortunately, 
one-dimensional models have an inherent weakness. The 
LSP tends to be highly dependent on laser beam geometry, 
radial convection, and radial diffusion transport of momen- 
tum and energy in the flow. Keefer and his co-workers'* 12 
developed a simplified semi-two-dimensional model that con- 
sidered the radial energy diffusion, but assumed that the 
laser beam and the streamlines were parallel. Later, Glumb 
and Krier 13 extended that model by using a converging laser 
beam as well as more realistic gas properties. However, their 
model still did not consider the momentum equations, but 
instead assumed a known velocity field. Merkle and his co- 
workers * 4 - 15 used a contemporary numerical method to at- 
tack the full two-dimensionai problem. They succeeded in 
the calculation of a low-temperature LSP (peak temperature 
about 4000 K) in a hydrogen flow seeded with cesium. It is 
feasible to extend their approach to a pure hydrogen plasma. 
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but the higher peak temperature (about 16,000 K) and much 
stronger nonlinear optical properties of the pure hydrogen 
and argon plasma can result in numerical stability problems. 

Recently, Jeng and Keefer 1 obtained a numerical solution 
for a rigorous two-dimensional model with realistic gas prop- 
erties for the LSP in an argon flow. Their results agree well 
with the experimental data 2,3 in the prediction of plasma 
position, size, shape, and peak temperature. Jeng and 
Keefer 4 have also studied the effect of laser power, flow 
rate, static pressure, and optical geometry on the hydrogen 
LSP. They found that the zero radial velocity assumption 
used by Glumb and Krier 13 is not adequate for this problem. 
Their work also indicated that there are no fatal errors in the 
laser propulsion concept and that the ethical parameters, 
e.g., plasma size, position, radiation heat loss, and energy 
conversion efficiency, can be controlled using an appropriate 
optical arrangement and a suitable flow configuration. 

The purposes of this work were to evaluate the full two- 
dimensional laser-sustained plasma model developed by Jeng 
and Keefer 1 against experimental results from Welle et al. 3 
to compare this model prediction with that from the constant 
axial mass flux, zero radial velocity, laser-sustained plasma 
model proposed by Glumb and Krier 13 (referred to in this 
paper as the “quasi-two-dimensional” model), and to ex- 
plore possible future model improvement. In the following 
section, the theoretical method used in the study is briefly 
described. This is followed by discussions of model calcula- 
tions and of comparisons between the model calculations 
and experimental results. Additional details concerning the 
theoretical method can be found in Ref. I. 


j i X (mm) * 




Index of Refraction » 1.491 
Max. Thickness = 39.9 mm 


Theoretical Model 

The object of this work was the verification of the full 
two-dimensional model, and the simple flow geometry as 
shown in Fig. 1 was adopted. The domain is axisymmetric, 
and a cylindrical coordinate system was employed in the 
model. The argon enters from the bottom with a uniform 
velocity and temperature distribution and an originally col- 
limated 10.6 /xm C0 2 annular laser beam is focused into the 
gas by an appropriate lens. 

It is assumed that the flow is laminar, low Mach number, 
steady state and axisymmetric and has variable thermo- 
physical and optical properties. The pressure of the flow is 
relatively high, so that local thermodynamic equilibrium 
(LTE) can be assumed with little error. Therefore, the 
plasma can be described by a single temperature and its in- 
trinsic properties are only a function of this temperature and 
pressure. Thermal radiation is divided into two parts: an op- 
tically thin portion in which all radiation escapes from the 
plasma and an optically thick portion that can be described 
as a diffusion process. Since the Mach number for the flow 
is small, the kinetic energy and viscous dissipation are 
neglected in the energy equation. Following these assump- 
tions, the equations of conservation of mass, momentum, 
and energy for the flow can be written as 

d{pu) 1 d(rpv) 

— — — + 0 ( 1 ) 

dx r dr 

8 (puii ) 1 8{rpuv) 8p 

dx r dr dx 

♦£( 4 ) 4 ±(-*)] 


Fig. 1 Sketch of the test configuration. 



10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 


axial distance (mm ) 

Fig. 2 Isothermal contour plot and mass flux vector plot of a 
calculated plasma (isothermal lines starting from 1000 K with 1000 K 
increments). 
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where I, is the local laser intensity described by 


(4) 


d(puv) 1 d(rpvv) _ dp 
dx + r dr dr 



(5) 


[ d / dv \ 1 d ( du \1 2pu 

L dx v dx / ' dr \ r dr ) J r 2 


( 3 ) and s, is the distance along the laser ray propagation direc- 

tion. The laser beam path through the optical system was 
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calculated by geometric ray tracing. The diffraction of the 
laser beam due to the finite aperture of the lens and the ef- 
fects due to the inhomogeneous refractive index within the 
plasma are neglected in the calculations. 

In Eqs. (2) and (3), the viscous term related to dilation has 
been omitted for computational efficiency. A few computa- 
tions were made with the dilation term included and the in- 
fluence of this term on the computed flowfield was found to 
be negligible. 

The following properties of argon were used in the 
analysis. Thermodynamic properties (enthalpy, specific heat 
at constant pressure, and density) were interpolated from the 
tabulated values given by Drellishak et a!. 16 The viscosity be- 
tween 5000 and 20,000 K was obtained from De Voto 17 *' 8 
and below 5000 K the viscosity was interpolated from the 
graph given by Vargaftik and Filippov. 19 The value of the 
intrinsic thermal conductivity for temperatures above 5000 K 
was adopted from De Veto 17 * 18 and for the low temperature 
region, it was obtained from Vargaftik. 20 The absorption 
coefficient at 10.6 pm was calculated using the expression 
given by Kemp and Lewis, 21 whose analysis included 
photoionization as well as electron-ion and electron-neutral 
inverse Brehmstrahlung using the Gaunt factors given by 
Karzas and Latter. 22 The expressions for radiation loss were 
obtained from Kozlov et al., 23 who considered both line and 
continuum emission in the optically thin limit as a function 
of temperature and density. 

The use of the diffusion approximation for the radiation 
transport is not trivial; it is not a gas property and it depends 
upon the actual size and detailed temperature field of the 
plasma. If, according to the study of Kopainsky , 24 the size 
of the argon plasma is greater than 5 mm and the pressure is 
higher than 3 atm, then the radiation contribution below 200 
nm represents a true optically thick case and the radiation- 
induced thermal conductivity no longer depends on the size 
of the plasma. The size and pressure of the plasmas under 
consideration are comparable to these limits. The theoretic- 
ally calculated value of radiation-induced thermal conductiv- 
ity 24 was used for static pressures greater than 3 atm, and an 
interpolated value between the experimental results 25 at 1 
atm and the theoretical calculations 24 at 3 and 5 atm were 
used for low-pressure (<3 atm) flows. 

Finite-Difference Solution 

The finite-difference solution procedure employed in this 
work is an adaptation of the method described in detail by 
Patankar. 26 The detailed description of the numerical pro- 
cedure used in the current calculations is discussed by Jeng 
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Fig. 3 Comparison of the calculated temperature contour between 
full two-dimensional and quasi-two-dimensionai model (contour 
lines starting from 10,500 K with 500 K increments). 


and Keefer 1 and only a brief description is presented in this 
section. The SIMPLE 26 - 27 algorithm, which uses a conser- 
vative form of finite-difference equation, primitive variables, 
and staggered grids, was adopted to solve the coupled mass, 
momentum, and energy conservation equations [Eqs. 
(1-4)]. 

At a temperature of 300 K, argon is transparent to the 
laser irradiation, and the flow will absorb no energy. If the 
numerical calculation starts with a 300 K initial temperature 
distribution, a flow with uniform 300 K temperature will 
result. This solution is trivial and no plasma will be 
predicted. To obtain a (stable) LSP, an initial high- 
temperature zone around the focal point was assigned 
(resembling a spark ignition) where the size, position, and 
temperature of the “spark” were chosen by trial and error. 

The computer code used in this study is similar to that of 
Jeng et al.'- 4 Generally, the convergence of the numerical 
solution was achieved in 2000 iterations and the CPU time 
requirements are approximately 40 s for each iteration on a 
Masscomp-500 computer. 

Results and Discussions 

In the calculations, the wall and gas inlet temperature were 
assumed to be constant at 300 K and the inlet velocity to be 
uniform. The CO z laser operated at a wavelength of 10.6 /xm 
and the collimated laser intensity distribution was that used 
by Jeng and Keefer.' It is an annular beam with peak inten- 
sity at a radius of 25.5 mm. The dimensions of the 8 in. 
focal length lens, which is one of two lenses used in the 
calculations, are shown in Fig. 1. Another lens of 12 in. 
focal length is similar to the 8 in. focal length lens, except 
for a maximum thickness of 39.9 mm and a radius of cur- 
vature of 143.2 mm. 

The fundamental physical processes occurring within the 
LSP have been discussed by Jeng and Keefer.' The validity 
of the present model was established by calculations for two 
typical LSP’s. The calculation included the interaction be- 
tween the laser beam and the plasma, thermal radiation 
transport, and fluid convection. In the present paper, the 
verification of the model has been extended to a wider range 
of gas inlet velocities, gas static pressures, and laser power 
using two different lenses of 8 and 12 in. focal lengths. The 
flow conditions were varied at a static pressure of 1. 3-4.0 
atm, the flow inlet velocity at 0.4-4.5 m/s, and the laser 
power at 261-967 W. 

Full Two-Dimensional and Quasi-Two-Dimensional 
Model Calculations 

The temperature gradient within an LSP is very large, 
resulting in very high-density gradients and thus a very com- 
plex flow in the plasma. The most sophisticated quasi-two- 
dimensional model assumes a zero radial velocity, leading to 
a constant axial mass flux assumption. This simplifies the 
coupled nonlinear momentum and mass conservation dif- 
ferential equations so that the energy equation is the only 
differential equation to be solved numerically and the stabil- 
ity of the numerical solution is easily achieved. From the 
results of this work, the CPU time required for the quasi- 
two-dimensional model is at least of an order of magnitude 
less than that required for the full two-dimensional model. 
Figure 2 shows a typical calculated plasma from the full two- 
dimensional model. The isotherms are shown in the top half 
of the figure and the mass flux vectors in the bottom half. 
The laser focal zone is located 50 mm from the inlet and the 
directions and the magnitudes of the local mass flux, at the 
vector origins, is proportional to the directions and lengths 
of the vectors. The detailed heat-transfer and momentum 
transfer mechanisms for this LSP have been discussed by 
Jeng and Keefer,' so only a brief discussion is presented 
here. The LSP is positioned near the laser focal zone, and 
the peak temperature is 16,390 K. In the high-temperature 
region, where the aixal velocities were much larger than 
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elsewhere, the axial mass flux was very small compared with 
the inlet mass flux. It was found that the plasma produced a 
local high-pressure zone near the leading edge of the plasma, 
which forced the flow to have outward radial velocities near 
the centerline. From this result, it is clear that the constant 
axial mass flux assumption in the quasi-two-dimensional 
model is inadequate. 

A comparison of the temperature distribution calculated 
using the full two-dimensional and the quasi-two-di- 
mensional models is shown in Fig. 3. The figure shows only 
a portion of the calculation domains near the laser focal 
zone and isotherms with temperature higher than 10,500 K. 
The innermost and outermost laser rays are also plotted in 
the figure. It was demonstrated in the earlier paper 1 that the 
calculations from the present model agreed well with the ex- 
perimental data with regard to LSP position, size, and 
temperature distribution. The quasi-two-dimensional model, 
compared with the present model, predicts a smaller size for 
the plasma and a higher peak temperature. Although the 
location of the predicted maximum temperature is the same 
for both models, the detailed temperature distributions are 
considerably different. The quasi-two-dimensional model 
calculated a very steep temperature gradient and a short 
high-temperature zone upstream of the laser focal zone and 
most of the plasma was positioned downstream of the focal 
point. From the mass flux distribution shown in Fig. 2, it is 
clear that the quasi-two-dimensional model overestimated the 
axial mass flux in the high-temperature zone. In this case, 
the calculated plasma using the quasi-two-dimensional model 
is similar to the predictions that would be obtained with the 
full two-dimensional model having increased inlet velocity. 

Model Verification 

The model has been compared with a wide variety of ex- 
perimental results given by Welle and his co-workers 2 - 3 in 
order to learn the capability and limitations of the model. 
Since the spectroscopic measurements were limited to 
temperatures greater than 10,000 K, the following figures 
contain only those isotherms within that temperature range, 
and a limited portion of the calculation domain is plotted in 
order to gain a clearer comparison. Due to space limitations, 
only a few typical comparisons are presented in this section. 

Inlet Velocity Effect 

The performance of the model for two different inlet 
velocities, using the 8 in. focal length lens, is demonstrated 
in Figs. 4 and 5, respectively. The static pressure and laser 
power of these two LSP’s differ by less than 1% and the in- 
let velocity of the LSP in Fig. 4 is about seven times that in 
Fig. 5, The model calculates the size and position of both 
plasmas well. As the inlet velocity increases, the radius of the 
plasma decreases, and the plasma is located further 
downstream. The local maximum temperatures within the 
low-velocity LSP are predicted by the model. One is on the 
laser path upstream of the focal zone and another is at the 
focal zone. Both calculated maximum temperatures are 
about 1000 K higher than the measured peak temperature 
(14,500 K). Only a single maximum temperature was 
measured, which was located neither on the laser path nor at 
the focal zone. The predicted local peak temperature at the 
focal zone, which was not observed in the experiment, may 
result from the neglect of the diffraction and refraction of 
the laser beam. Although it has been demonstrated ex- 
perimentally 2 - 3 that this assumption leads to only a small ef- 
fect on the measured global laser power absorption, the 
predictive model can suffer more serious effects. If the 
refraction of the laser beam within the LSP is considered, it 
results in a larger focal spot size and reduced laser intensity. 
The spot size for the refracted beam is highly dependent on 
the L3r> position and size. For a low-velocity LSP, most of 
the plasma is located upstream of the focal zone, making the 
focal spot larger. When diffraction is also considered, the 



axial distance (mm) 


Fig. 4 Temperature contour plot of a high-velocity plasma (contour 
lines starting from 10,500 K with 500 K increments). 



axial distance (mm) 


Fig. 5 Temperature contour plot of a low-velocity plasma (contour 
lines starting from 10,500 K with 500 K increments). 



Fig. 6 Temperature contour plot of a low-velocity plasma with the 
12 in. focal length lens (contour lines starting from 10,500 K with 
500 K increments). 
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axial distance (mm) 

Fig. 7 Temperature contour plot of a plasma with high laser power 
(contour lines starting from 10,500 K with 500 K increments). 



axial distance (mm) 


Fig. 9 Temperature contour plot of a low-pressure (1.5 atm) 
plasma (contour lines starting from 10,500 K with 500 K 
increments). 



axial distance (mm) 

Fig. 8 Temperature contour plot of a plasma with low laser power. 


diffraction limited focal spot is approximately twice the size 
of that calculated from a geometric ray trace. The combina- 
tion of laser beam refraction and diffraction and the longer 
absorption pathlength ahead of the focal zone result in an 
actual laser intensity that is smaller near the laser focal zone 
than that which is predicted by the model. Thus, the local 
temperature maximum will not be generated there in the low- 
velocity LSP. This may explain why the calculations are not 
consistent with the experiments in this detail. 

Since the size and temperature of these LSP were well 
predicted, it is not surprising that the calculated thermal 
radiation power loss (316 and 445 W for high- and low- 
velocity, respectively) agreed with the experimental results 
(326 and 392 W for high- and low-velocity, respectively). The 
model prediction for the transmitted laser power in the high- 
velocity LSP also agrees reasonably well (110 and 159 W for 
prediction and experiment, respectively). However, the 
predicted transmitted laser power (91 W) in the low-velocity 
LSP is much smaller than the measurement (246 W). This is 
because the predicted second local maximum temperature 
near the focal zone absorbs an unrealistically large amount 
of additional laser power. 

Lens Effect 

Isotherms are shown in Fig. 6 for an LSP with operating 
conditions similar to those of Fig. 5 except for the use of the 



axial distance (mm) 

Fig. 10 Temperature contour plot of a high-pressure (3.0 atm) 
plasma (contour lines starting from 10,500 K with 500 K 
increments). 


12 in. focal length lens instead of the 8 in. lens. Both calcula- 
tion and experiment show that the plasma moves further 
upstream than in Fig. 5 and the model predicted the 
upstream position and the radius of the plasma well. Again, 
the model predicts a local maximum temperature near the 
focal zone, which was not measured in the experiment. The 
effect of laser beam refraction and diffraction, which was 
discussed in the last section, is probably responsible. Similar 
to the 8 in. LSP, the thermal radiation from the LSP is well 
predicted and the transmitted laser power is underestimated. 

Laser Power Effect 

Figures 7 and 8 show the effect of laser power on the 
structure of the plasma. Both cases use the 12 in. focal 
length lens with similar flow conditions. The calculations are 
quite consistent with the experiments. As laser power in- 
creases, the plasma becomes larger and moves upstream. For 
the 10,500 K isotherms, the model underpredicted the length. 
However, the experimental results exhibit more noise in these 
isotherms than in the higher-temperature isotherms; if the 
11,500 K isotherms are used, then the model predicts the 
length of this isotherm within 10%. The calculated positions 
of both plasmas are about 1 mm downstream of the ex- 
periments. After carefully reviewing the experimental data, it 
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Fig. II Temperature contour plot of a low-pressure (1.5 atm) 
plasma using 50% of radiation-induced thermal conductivity (con- 
tour lines starting from 10,500 K with 500 K increments). 


was found that the on-axis maximum temperatures are not 
on the laser beam path (using real ray tracing and neglecting 
the inhomogeneous refractive index effect). If the experimen- 
tal results are shifted downstream by 1 mm, which is within 
the accuracy claimed by Welle et al., 3 the peak temperatures 
are located on the laser path. In this case, the model predic- 
tions agree well with the experimental LSP positions. The 
predicted global properties are also in reasonable agreement 
with the experiments: the transmitted laser power is within 
10%, and the thermal radiation from the LSP is under- 
predicted within 30%. 


Static Pressure Effect 

All the LSP’s shown in the previous sections have either 
2.0 or 2.5 atm static pressure and the model performed well 
at these pressures. In this section, the effects of the static 
pressure are presented. The intrinsic properties of plasmas 
are highly pressure dependent, especially the optical proper- 
ties. The absorption coefficient of the plasma at the laser 
wavelength and the thermal radiation from the plasma are 
approximately proportional to Figures 9 and 10 show 
the plasmas at 3.0 and 1.5 atm, respectively. Both plasmas 
have similar laser power levels and inlet mass flow rates (but 
different inlet velocities). In the low-pressure case, the 
predicted plasma size is much smaller than the experiment. 
In the 3 atm case, the model provides an excellent 
prediction. 

After reviewing the model assumptions, we concluded that 
the treatment of the optically thick thermal radiation was the 
weakest link and was probably responsible for the failure of 
the model in the low-pressure case. As mentioned before, the 
validity of the diffusion approximation for short-range 
radiation is dependent on pressure and the size of the plasma. 
As the pressure decreases, the required size of the plasma in- 
creases. For the present calculation, the size of the 1.5 atm 
plasma is below the minimum size required for the valid use 
of the diffusion approximation found by Kopainsky. 24 
Although the diffusion approximation for short-range ther- 
mal radiation is not valid for the low-pressure plasma, a sen- 
sitivity study using this approximation was performed to in- 
vestigate the sensitivity of the solutions to radiative 
transport. For small, low-pressure plasmas, the radiation- 
induced thermal conductivity should be smaller than that for 
a large, low-pressure plasma. 

A calculated result, where the radiation-induced thermal 
conductivity was reduced to 50% of the value used for the 
calculation in Fig. 9. is illustrated in Fig. 11. The size of this 


plasma is much larger than the previous calculation and 
agrees well with the experimental results. The role of a 
change in thermal conductivity on the results can be 
understood qualitatively. As conductivity decreases, the con- 
duction heat-transfer loss from the portion of the plasma 
along the laser path is reduced, which increases the plasma 
temperature along the laser path. Due to this increased 
temperature, the absorption coefficient at the laser 
wavelength increases and the plasma absorbs more laser 
power, which in turn generates a larger plasma. Because the 
exact solution of radiation heat transport will not be easy to 
achieve in the near future, the optically thin and thick ap- 
proximations for the radiative transport process should be 
retained for the near term. Although in some pressure and 
dimension ranges the optically thick approximation is not 
quite valid, this concept can still be utilized with some ad- 
justment for the magnitude of the radiation-induced thermal 
conductivity. 


Conclusions 

A full two-dimensional model for laser-sustained plasmas 
has been numerically solved for a wide range of argon flows. 
The calculations were compared with experimental tempera- 
ture distributions and it was found that the model performed 
well at static pressures greater than 2 atm. The plasma posi- 
tion and size are well predicted and the behavior of the 
predicted plasma structure using different inlet velocity, 
focusing lenses, and laser power agreed well with experimen- 
tal results. 3 

At lower pressures, the model underpredicts the plasma 
size. It was recognized that the diffusion approximation for 
short-range thermal radiation is probably responsible. When 
50% of the optically thick, radiation-induced thermal con- 
ductivity at a pressure of 1.5 atm is used, the calculations are 
in reasonable agreement with experiments. 

The global thermal radiation heat loss from the LSP and 
transmitted laser power are well calculated except for two 
LSP’s. In these two LSP’s, two local maximum temperatures 
were predicted. However, the experiment did not show the 
same two maxima. We concluded that neglecting the laser 
beam refraction and diffraction is responsible and recom- 
mend that the model should include these effects in future 
studies. 

A comparison of calculations from the present model and 
those from a constant axial mass flux model was also per- 
formed. The axial mass flux field calculated from the present 
model differs considerably from the constant axial mass flux 
assumption in the simplified model. The characteristics of 
the temperature fields calculated with these two models are 
different and the results of the simplified model are similar 
to the calculated results from the present model using a 
much higher inlet velocity. 

Acknowledgment 

This work was partially supported by U.S. Air Force Of- 
fice of Scientific Research Grant A FOSR-8 3-0043, the Pro- 
gram Manager was Dr. Robert Vondra. 


References 

'Jeng, S-M. and Keefer, D. R., '‘Theoretical Investigation of 
Laser-Sustained Argon Plasmas,” Journal of Applied Physics , Vol. 
60, Oct. 1986, pp. 2272-2279. 

2 Keefer, D. R., Welle, R. P., and Peters, C. E., "Power Absorp- 
tion in Laser-Sustained Argon Plasmas," AIAA Journal , Vol. 24, 
Oct. 1986, pp. 1663-1669. 

3 Welle, R. P., Keefer, D. R., and Peters, C. E., "Laser-Sustained 
Plasmas in Forced Convective Argon Flow, Part 1: Experimental 
Studies," AIAA Journal , Vol. 25, August 1987, pp. 1093-1099. 

4 Jeng, S-M. and Keefer, D. R., "Numerical Study of Laser- 
Sustained Hydrogen Plasmas in a Forced Convective Flow," AIAA 
Paper 86-1524, June 1986. 


ORIGINAL PAXZ is 
OF POOR QUALITY 



1230 


JENG, KEEFER, WELLE, AND PETERS 


AIAA JOURNAL 


5 Mcrklc, C. L., “Laser Radiation to Supply Energy for Propul- 
sion,” AIAA Progress in Astronautics and Aeronautics: Orbit - 
Raising and Maneuvering Propulsion: Research Status and Needs , 
Vol. 89, edited by L. H. Caveny, AIAA. NY, 1984, pp. 73-94. 

6 Keefer, D. R., Elkins, R. f Peters, C., and Jones, L., “Laser 
Thermal Propulsion,” AIAA Progress in Astronautics and 
Aeronautics: Orbit-Raising and Maneuvering Propulsion: Research 
Status and Needs , Vol. 89, edited by L. H. Caveny, AIAA, NY, 
1984, pp. 129-148. 

7 Raizer, Y. P., “Subsonic Propagation of a Light Spark and 
Threshold Conditions for the Maintenance of Plasma by Radia- 
tion,” Soviet Physics— JETP, Vol. 31, 1970, pp. 1148-1154. 

8 Jackson, J. P. and Nielsen, P. E., “Role of Radiative Transport 
in the Propagation of Laser-Supported Combustion Waves,” AIAA 
Journal . Vol. 12, 1974, pp. 1498-1501. 

9 Kemp, N. H. and Root, R. G., “Analytical Study of Laser- 
Supported Combustion Waves in Hydrogen,” Journal of Energy , 
Vol. 3, 1979, pp. 40-49. 

10 Keefer, D. R., Peters, C.E., and Crowder, H. L., “A Re- 
examination of the Laser Supported Combustion Wave,” AIAA 
Journal . Vol. 23, 1985, pp. 1208-1212. 

u Batteh, J. H. and Keefer, D. R., “Two-Dimensional Generaliza- 
tion of Raizer’s Analysis for the Subsonic Propagation of Laser 
Sparks,” IEEE Transactions on Plasma Science, Vol. PS-2, 1974, 
pp. 122-129. 

l2 Keefer, D. R., Crowder, H. L., and Elkins, R., “A Two- 
Dimensional Model of the Hydrogen Plasma for a Laser Powered 
Rocket,” AIAA Paper 82-0404, Jan. 1982. 

,3 Glumb, R. J. and Krier, H., “A Two-Dimensional Model of 
Laser-Sustained Plasmas in Axisymmetric Flowfields,” AIAA Paper 
85-1553, July 1985. 

u Merkle, C. L., Molvik, G. A., and Choi, Y.-H., “A Two- 
Dimensional Analysis of Laser Heat Addition in a Constant Absorp- 
tivity Gas,” AIAA Journal . Vol. 23, 1985, pp. 1053-1060. 

13 Mcrkle, C. L., Molvik, G. A., and Shaw, E. J.-H., “Numerical 
Solution of Strong Radiation Gasdynamic Interactions in a 


Hydrogen-Seedant Mixture,” AIAA Paper 85-1554, July 1985. 

l6 Drellishak, K. S., Aeschliman, D. P., and Cambel, A. B., 
“Tables of Thermodynamics Properties of Argon, Nitrogen, and 
Oxygen Plasmas,” AEDC-TR-64-12, 1964. 

|7 DeVoto, R. S., “The Transport Properties of a Partially Ionized 
Monatomic Gas,” Ph.D. Thesis, Stanford University, Stanford, 
CA, 1965. 

18 DeVoto, R. S., “Transport Properties of Ionized Monatomic 
Gases,” The Physics of Fluids, Vol. 165, 1973, pp. 615-623. 

I9 Vargaftik, N. B. and Filippov, L. P., Thermal Conductivity of 
Gases and Liquids ( Data Book), Standards Press, Cincinnati, 1970. 

^Vargaftik, N. B. t Tables on the Thermophysical Properties of 
Liquids and Gases , Hemisphere Publishing, Washington, DC, 1975. 

21 Kemp, N. H. and Lewis, P. F., “Laser-Heated Thruster Interim 
Report,” NASA CR-161665, 1980. 

^Karzas, W. J. and Latter, R„ “Electron Radiative Transitions 
in a Coulomb Field,” Astrophysical Journal Supplement Series , 
Supp. 55, Vol. VI, May 1961. 

23 Kozlov, G. I., Kuznetsov, V. A., and Masyukov, V. A., 
“Radiative Losses by Argon Plasma in the Emissive Model of a 
Continuous Optical Discharge,” Soviet Physics— JETP, Vol. 39, 
1974, pp. 463-468. 

24 Kopainsky, J., “Strahlungstransportmechanismus und Trans- 
portkoeffizienten im Ar-Hochdruckbogen,” Zeitschrift fur Physik , 
Vol. 248, 1971, pp. 417-432. 

25 Bues, J., Patt, H. J., Richter, J., “Uber die Elektrische Leit- 
fahigkeit und die Warmeleitfahigkeit des Argons bei Hohen 
Temperaturen,” Zeitschrift fur Angewandte Physik, Vol. 22, 1967,- 
p. 345. 

^Patankar, S. V., Numerical Heat Transfer and Fluid Flow, 
Hemisphere Publishing, Washington, DC, 1980. 

27 Gosman, A. D., Pun, W. M„ Ruchal, A. K., Spalding, D. B., 
and Wolfshstein, R., Heat and Mass Transfer in Recirculating 
Flows , Academic Press, London, 1969. 


B-39 



APPENDIX C 


Computer Code 



ORIGIFif.L Ff, 
OF, POOR QO 


« « -p 

* * «-i G 

* * «M | fit 

* * I I * 

* # H H CN 

« « ZZl* • 

* * I I o I 

« * gnvo 


s : S 5 oi 

« « zzfi! 

« « 

* * 

* * 

« « 

* « 

« « 


i?S?a 


x' u y 

+ fa w 
^-»in >< o 
fa • t K 
£7 * 

Se? h> I 

• fcH S -r"s 

8fga 

• (D H H 
£-CT> Z2 

wh^ x rH 

£*AAA 

fa fa 2 03 2 
M fa < CQ < 


ffl cci* ^ 

CQ CQ h) - 


* * 

£ 

I M 

h> h> M 2 

?7>? 

M W K *-s 
221 ^ 


fTlH |g? 

- -W | 

H H 2 

S 44 » 
H H !C ^ ^ **> » '-' 

t 5 > + ^ ■ — >VH w <•*. 

D 2 Z fa 

+ + b 25 

v v /s A, fa A 

t? l?HZ> * • * • 

grfs???I 
221 £ ^ 1 . *1 
5 *>l 

«ii.{?BM?SaB. 

o'I^mdS.SSE 1 

« m m 2 wd cu a a i 
| M fa&TSTfa 

» D > H MMWM; 


\ 

CQ hi 

g? H 

hi I 

►3 m w 

7 | f 

wt X 

Z ^ *-) 

tsti 

3‘i® 

2 r*— - <-si 

z ?B £ 

* + H *«*H 

'flfl 

S888 


fa g fa 

£ P I w 
fa ■* ^-w 

O' (3 H fl< 

* P fa 
l *— • fa fa 


fa l "3 hi hi 

M tH M 

h as z z 

M fa fa fa 
Sfafafa 
fa H HI M 


2 O fa CO 

fa»fg 

H <s^O # 

fa hi hi to <* 

rsr s '5a 

2 ^ l 

^-Sfa « A 

U fa fa K M 
fa W M * fa 


J~a| - 

Cu °ix • 2^0 ^-. 

•H H M - ~© r-s 

Se§ - 8 ^ 8 ^ 

•PE^hhmh 
Pb 2 « 

^glpA^gb 

20 2 *j 8 



a 

M 

>» 

»J 

W 

« 

so 

« 


* CQ CQ H 2 5 
rH CQ CQ WM > 

£®£>Sl 

O fa fa fa Q < 

u WHHOI 


f r' 

a g 

CQ - 

8 ? | 

Bsc* 

B « *Ir£ 

ill 


agR\ 

sis saga 

El| -gig 

- - p Ph *• *- 

2 o> fa h P P 

38 *. i?a s .i 5 

ESS agEe: 

- ^fa fa W — - M 

<— I (N OJ - > ^ 

S 882 g°.e s .sa 

) fa fa fa fa P3 >*% «. 

> fa fa 03 fa S ^ P 

* " ^ fa ^ M s?3 

I'E S °i-S M SH g« 

.egsgEi^ra 

: a 

H(Nn^r\^ •*. ^ 
lUUBHfeE^ 

lilsSbiBi^^ 


-IH B 

fa m w 

h ~ W 

e* a 

W /% \ 

as P r 

es P 

P fa w 

M -.^>2 

ffjS 

x .sa? 

S>iS 

*.a 5 v 

RB|g 

w PpS 
2^P5^ 


‘Ilia! 


:|IIIl 

!8s9B» 

p.M'aas 

g^a*' 


ills 


SVgDg 


S- fr. 

H W fa W 


B c g-gePBS S 

aS a g5gaa|i I 
&iaps^8?aa a .c |g 

!*e^°eg=o'--p l i is 

ig»SPx E ;y5gg f iS la 

I H ^ ^h < H Ch 2 

^ggE;S<§SSSa||ig 

2 2 Eh M <«H£ 

MObp^< \ '•ooh <<P5 


12 2 2 2 2 ! 


!Sft w||! 


Eh m -fa 2 

zzrh ^E^sso > 

M o b C* \ -CD <T« rH 

-u -\x 2 o?° Q Q Q 

\P qs< 3 qsoqqq 

^ HHg \\\ 


gZH -P g g 

ill lilt 

05 OS - -fa 05 W 

ggS.zSSgS 

§^§ei|.s§. 

Z jS Z 2 Z 


«Hrg*n ^*inso\o rH 


uuu 


C-l 




£ uuu b o o u uuuuu u u u uu 



C-2 




CN *fc 

Ml 

U ml i> r** <x> r- en 

<$ tn in in ^ o 

S V v vo oov o+ t- v uun OUUUH 


C-3 



^ 8<* 


I 

cl 


CO 

00 

OS 


r>F POOR QUALITY; 



S 


S2 


ID 

(N 


I 


£ 


S 2 < 


8 

r 

H H^) hj 
ZZH M 

$8? ¥ 

/N ^ /«»S 

HHbhH M 

S £ - -+ + in 

* y5 HH ^ in ^ * 

I - ^Z Z m * m 

■ — * M M H KJ w,- 

^,. I I SSqs ^bc *h 

M 5 ***+*1-5 

^z«(Ninmzb2 ' 

* MOO • • z - Z M 

- PP,H ?| 

« — w 


+ 


in 

03 


* 

?3 

3 8 


S3 


£"> »-5 ~ 


8 

« 

i-> 


M £ 

> « 

iH #-n 

I I 

<-n. - 

►i M 

M 06 

«■— <^* 
OS M ^ 


+ 

•-3 


z 


*7) m 3E <-* 

i i S3 

K 5? iS 

51 
§| 



^ H«M -+ in SO 


! + i 


H b H 2 b 

g£g£gd 

• 2 • Z • Z 

S8&8£8 

6. T h T ta T 


g g i ? 3 * 


Si £ 

£? ¥ 

OS ^ ^ 

m ^ -in i 

~ J iH W • M 

)+w*w 

: t -3 os os 

O M ** *^ 

“zing 

sggig 

hTT — 7 

HINH WN 

1353 

03 M 


OS 

O 

in 



M in 

<d 00 

£ <J O m 


so 

00 

in OO O 


c-k 



« 

« 

« 

* 

« 

« 

« 

* 

« 

* 

« 

* 


CO 

03 

o> 


Z 

03 

« 


*3 

H 

+ 


X S* 


rH 

+ 


*"l 

«H 

+ 

M 

I 


*1 

H 

+ 


8 


►0 


0 o 

1 I 


*n ►o 


• Z Z * 

•teat • 


»-j rr 9 9 


►o S* S Mgrsw^SlSw 

£ f B^dd^a 

a: KH ^H^hh »SHhhH 

* * z z z z 

#-*,*-**/-*<«*« «* M « »R 

H^H^ZHZH I I H Z QQZH 

gds£g|s|55ss^^sg 

EMiUtibabiSSbb tb b* 

HHHP3HHH05<<<HH M M 


0 - 1 * £ *-*o-s 



S* -fee 

H z h V os 


+ S 
? • 
l * 


B 3 S 

J-*§ A u if» 


Cu Cu Hinomo 

+ + S • + • + 

o-* • + O 0 -*O o-N 

H>rq ^ 

A b o-*rH /— >"3 

H *“'» «*. 0 -N | O-S *. 

as mir 


-I 0* OS H 

+ W'CU 

^ \ tL 

- 0-* /~* + 

M Hj P** ^ 


§ a § a 



►o 

rH 

I 


05 05 ! 

I I I 

M 

a> O' 

-# in 

ro ro 

28 


A 


52*1 

!?v»$ 

0 * D H 

03 0) OS 


4S 

M (-J 

ST m 

03 w 

*5 

M -* 

r 

^z © 

h h n 

+ * rH 
»"3 * — «■ 

SBh‘ 

?Sa 


SSi* 


77 

M l-J 


; f 3 /non 

'4b 


sr sr 


88 


03 

0 - 0 * O 

u u ^ 


C-5 





************ CENTERLINE BOUNDARY CONDITIONS ***************** SUBROUTINE CALCU 


& 


00 

00 

o 


3 


b 

pH 

I 


b 


H f 


D p - * 

I b 




si 

£V 

3 + 

W M 

i& 

Eft? 

41 CN «*■ 

in . h 

CN « | 




D ^ 

J-V 

+ ^ 

32 

|? 

b M 
m D 

D * —s 
'-in pH 
* 09 I 
in • b 
cn * - 


+ pH pH I 

pH^pH 
I pH pH + 
H I I M 
wM M w 

? DO? 
— » | | ^ 
pH ^>p>rH 
+ pH H + 

b I + b 
J“b b - 
pH ^ '-pH 
+ H H + 
M + | M 

2 r w M ^ 

SssS 

« '— — 4t 

in * •♦« in 
• in in * 

* • ♦ * 
-N* * — ' 

W | ^ M 

*— b pH *•— 
H ' | pH 


* t* r** ^ 1 


SB 2 

CN CN 

M A 


88 


f? 4 ? 


riri^ 


iSf 


IV S' 

M M | 


I I M 


S??g 

rj p*"s M 

Z CN CN 2 

SSgg 

^ 3 bHH 

Du Dm Du Cm 


b 

pH 

I 


b 

M 

DO 


I 

b 


I U U 


ORIGINAL „ 

OF Pone *s 

O0ft QUALITY 


b 


* H U w O 
4 c w « W « 

« 2 in ea in 

* 8 o 8 o 

* * + « + 

* in p> 

: ???$ 

* p's b p-'pH 

4 * ~ ^p^ I 

* hHbH 

: ifig 


u 

tH 


p"»pH p-' b 




/"S 

s;u 

H + H 

2 ~ 5 T 

*5 Cb I 


\ 

b 


\ 

b 


I 

b 


b 

v 


£ 


t 


\ 

«p» 

b 

M 

oT 

« 

d 

pH 

►o 




2 z oa Da ^Ea o> i b 

UUUUbW i ~ ^ '- 

U U U CJ '>4 ^H pH Cb 

* * « * h * h l I « 

in in in in ^ m 

• * * -WbH ' w <- 

0 0 0 0 2 - I H CU w 

wW M M W * * 

gj ^ J ^ w JH? 

V . 41 p^<-»4« M HI p-* 


Cb — p- ' 


N 4 + 
s— — * 

* 1-3 

PS s 

■* ps b M 
M D 


•s a 


f?ffi 3 £a 3 slati 8 l? 

vis 

i h T 

aaalSSilg "i 'i H s 


^ 3^22 

M V ^TThTfcJi^TfelwffliM V 5 S 


fflDCQQQjgg 


CQW + + Wi 
‘ CN 


pH O 
O O 


3 a 

oT d 
w - 

1 M 

2 2 0 
< h cn 
+ l 


O 

pH 

+ 

W 

D 

4 t 

p-s 




?a*? i. 

H ^ 

l£s ni 

p-s HI W P^ V + 

h U 5 Qj Hj H PS 

> flj < >■ sppH 

ag|aa^ 3 ? 

ZZHPSPSHH^ 

oh 4 i rr n ^2 

as o o v w s H i 
OnnHhtuHw O 

§saa MH gsi 


uuu 


u u 


& 

? 

Ob 


00 

00 

o> 


O b 

o ^ 
in z 

43 


-S 6 ! 

VoS2i 

kill 

j - L O Cb Cb l 

> O icT rH CN 
H M O OS OS 
JWOOHOp 

5 § Jj i Cb a 


P i ^ 

*t*is 
2 ^ 0 * 
M « 2 

H S d 
Du Cb a. 

" S' Cb 
pH CN 


K 


e 

B 

S 

8 ? 


K 

-fe' | 

BBp* 

s^ph pn 

SIb^ 

ft*lsR 

O p^*> •— 

S E^p d 

jiaBi? 

a w ww 
4 CN ^ > '- 



2 g°.p Q .a ' 

> a- cb d a Eo d S S* \ 

tn dS S °^S B 2 H Eh h 
QS K O O H O s-W '-'W 
OpESuCbfeCbWpH'— 0 

ESWW^BS 8 

H{Nfnf\^ ' ' 5 

aoooHHM&aou^^p-^^' 

sssJiiiiS^? 

????igas? 


ft !3 8*5 

p Cb "S M 

8 M V 2 S 



2 2 2 W 3 S D 

s " 6 ' s ;SiP. 




pH pH pH pH pH CN m <*• in ID 


h 

» 


U O O 


C-t 



ORtG^AL 


OF POOR QOAUS'J 


BIJS 

fill! 


b* E- 

h h Cm h 


\UM 

:P.y 

:es5s 


* F M ' 

^|| |||| 
IiLilil 

m o Z'Z'Z* 

i sl^ili 


ss|g: 

§saa; 

^gg 1 

ss 


+ HH I 


tfSTt 

i-H «■% iH 
+ HH + 

w- 


t>>t 


-■ 

1 Hj -'-**-* »"J 

vr-4 v 

* M I ^ M 

* '-•h h ^ 

* r- 1 »■ I »H 

'^Z MM u 

I H Q H H Q 

I^ASSA 

<m§??S 


rSSSS 


:> hfabfa 

| | W HI M M 


'Ig’tg 

HU H U 
— w* 

2 in h m 

BdBcJ 
* + * + 
in ^tn ^ 

l “3 ^rH 
I-J H 1-3 M 


Z.'TZ’l 

^ ^ ^ H 


Z Z ta H ^2 

B888^> 

« * * * M « 

in in in in ^ ^ 


X X X X M r- 

mifj 

*-s ,*-* ✓—* h 
^ ^ ^ ^ 


mm 

88 


in 2 05 o 2 Q | 
hvo u$ 
oo ^ oh! 


m ^*h 

as A § § § § cu cu cs cu P 
rH - i Hj a, a. a, a. ^ o* - _ 
Z Z < ' - - - H(S(£^ E 
O O w © *— I w fn ^ os OS - os h 
MH O OS OS OS OS p p HO' 
wwfflHppppEEaslc 

feiiaissss^^ 

HZSShi H^n^^ 

3 H H U —00 OOQOOOUUUUli 

u q q os * * * * 2222s 

- - OSQSOSQSOS<<<<D 

H.. - 


po po^HZ b i4 hS hS < a « & « z z z z 
faOSQOHQFSPtHHiMaUUUUgSgg 

Maoaagb||| 8888 


IcKKKl 

! £ h eh v £ c 


m 5 a a : 
w 5 < w < 


i ^ 2 I 

EEfc^S’ 

os os os iii v ; 

555tfH. 

e? > a -5 ! 
agaezf 


222| Z* 

88l|i‘ 

Ed ta U W D 


o 

HO O O 

in o o -«r o 

tn m vo 


C-8 



ORIGINAL PAGE 13 

OF POOR QUALITY 


(0 

CL 


09 

09 

<5% 


M fa RWfH ^ 

fa CL 0. m H 

^ -a o, oi 

H <N 03 w 

a a or>o g 

li^ip^pl 

£ £P £ b PES P 

^PpaPpSp* 

iipi^s^s 

fa fa fa fa fa 05 rH l w 

^^§.3^85 


3 z S D fa aT 

^ 05 Q < 05 




•« in VO VO H H 03 


M 

S 


& 

fa 


M 

S 


* 


I 

r> 

rH 

I 

M 

f 

►9 

H 

+ 

M 

fa 


5? f 


i 

►9 


r> 

\ 


\ 


4 


t 


fa w w *• *-* 

Z * fa '-►9 

M ~ * * <- 

< >“) « 

2 - '*&**£' 

K> 


ht 

E oj ~Z j 


o' o a 'O 

rH « H M S. 


TsH 


e 

CQ 


Tp 


CN 

« 

* 




8 

* 


fa 

05 


OJ 


m 


S^gi 
_.i,SS* 
n 5 w ST 


Pft 


5T^ 


V> 

Sf 5 s 

ASii 

O - • 
" rH 
HH - | 


+ Sh 

*-*M — 

i«a 

M 

Sod 


►9 

rH 

+ 

W 

F 

*•*» 

►9 

v. 

W 1-5 

as 

+ ED 
•'"* 05 
H + 

I <"* 
►9 ^9 




I 


gg?aa 

ZZH ^ /*> 

x - h 5 

> na c* Z - - 

♦I I < M H 

> H >5 | «- ^ 

?aa 


a= ,^Mg| 


I + H* 
rH 03 


©H , 

O O — 

nrOHhhi 


88 

H O 

o o 


!88a 


;:aa~ 

r*s -+ I fa 
h5 M r./*«0 
>-H 1-3 Cfi 


ga*** 

S ^HH 0 ) 

>^>>a 

b w« « + 
HHbbBS 

8f ss| 

^*9 < < £ 

i ~rH 05 ^ 

> > W &3 
Q Q fa fa 


>< 

3 ^ ^ ^ ^ C 

H W H H Wfa 

j»| 

LJ fa ^•rn^.^rH fa « 

gPbbbb —5 
5S ^ ^+# b 

M H H H *-**-» *« 

: paaaapp* 

* — i i i i 

■fc fa «rx*— s «■■'* <r>. ED ^ Q 

* < *-3 >"5 >9 >73 05 D I 

« || ' v v v | | r-*.r, 

* ^ •"> -2 

* >• > !> " - M 

* wzwhsmhhw- 

5 a ' c ~~'’ 

* 


O fa 

o ** 
in z 



w 8 8 


£g 

88 

H O 
O O 
m m 



*» Jr 

O «OHMn 

HH O fa fa fa 

05 05 00 H O Q O 
~ * J fa fa fa 

--.C3=3 , 5^' a,a,ai 

►3 S s < m 

<HHU wGO 00 05 

U Q D K * * * 

- - fa fa fa fa 

Elf 


U V V 


C-9 



r' ' ^ "F * q r~ \ n T-- 5 

; •; • * i* •; * f . » 

V-, 1 I i-m *'■..+■ -J t!%-' 

OF POOR QUALITY 


n 

£ 

M 

§ 


O 

e« 


I As. P 


00 

00 

ON 


s 



< * < o ^ 

2 r*± g h> ^ 


a 6 


fiffi 


h wQ F a 


M I OD I 
] M M | I 

| /«» /N 


• ww 
[03OHZG0 ZZQO ww 
m m r i *h «» 'Haa, 
z in m m w w 

SS838d>Sm& 



M 

b 


Q O 3 

cl| 

is* 


& 

r> 

a 

R 


- ( M O 
h W W 00 H 


baa 
"» "• a a 

H CN 04 

| o 

BBSSaacga^o 
- - » oT a ^ 
aSSaggfHg 1-1 

SSSSapuaarH 
a a a a \\\\p 
C3 ?!! 3! >> 

H 



U U U 




00 

00 


S 


i 

n 


M «-» ***««» * 

w «H »M • 

X +HH I * 

p hi lb * 

I - hi hi " « 

*-» M «*• **rH ft 

rH I HH + « 

I M | | M « 

Pj wM Mw « 

=_ !ffE I 

+ 0 
or 

M M 

II 


- «-^ + * 


+ rld + 

1 « + * 

+ H H + 
M + | M 
wM M w 

mi 

ft W w« 

in « * in 
• in in • 


ro 

hi M 


& E~rE£ 


rOEn 

m m —in i-j 
o« l * 04 ^ 
• m in • »-• 
* w(N * I 


*-> + « 

*» <— I ►o *» 

M| M 
wl-3 M w 

• I «H 


. S * hi w 

i ^ - r 

3 ^h3 M H 
*«H - I 1 


M | M M ^ 

s2asS 


E 

8 


A*0A,i,g??g £ 


sasg 

oM 
04 04 Z 

gggs 

>3^HH 

b b b h 



b 

M 

I 


w> <>s w 

0308 

M U M U 


* O O O O 

* * + * *f 

* in 

: ?f?o 

* oh ^*h 

* ^ l 

ft h> M 1-3 M 

: i%i$ 


z '--H 

b ^ i «-* *»• 

M hi hi *“i «-t 


I ^ 

M M « M *C 

Z Z M H 2 

P^P 

in m in in ^ 

O O* O O M 


« # « — s# « « — * 

^*“3 

rs 

M *> ^ •*• *>. ^ 

H + H I 3 H »• M *H *-> 

w M W >V — M W I )-3 

tsr^r^sr^O 

►oS^wmAmw* D 

v*v> + *i>< og w 

+ n -cn * * < j 

m r> m * m m m « «-«.go y o 

* 

C+ Ci ^C kwmmw^ ^ 

+ M+^+M+l|iCl-3© 1-3 

oa wrq m cfl woj hi m n *-0 
« D « I « D * ^ w * M 04 M 

•Hw*M*w#M& 3 ^w w 

^«.^x^*.Zx8rE8 gS g 

OTOrOTOrioSo f? s. 

M DO M • M Z M • « Mp Q 1-3 

w< w^» W< W^* N wQ h3 1-3 1-3 V 

DN©I©NDI*Zm^^-^ m 

”3t~t3A~AS3S820«”tl S' 


!???!. 

h) hi hi rs 1 1 

m m m m a a 
w ww^u y 

3333 


* M 


0? 


D 5 

C/5 + 
Ofl 


oso N rgr N Miigs5i 

rsr«o|o^iv*??s t 

M • M ■ M ■ M . y ft | W 

wO«o«O w Ou h h d 

d a d fl © a d a wh, ^ a - vp 

w . w . w . w . I - -MMM-.mI 

1 M J M | ^ I W M M w ww wA g g Z \ £ 

3 &T 3 b 3 b 2 b B z a&r ¥ 2 ow S B S 5 s j 

Nh^hNhNhhUWhJOWhHZZZS 3 

2 V 8 f 888 s 


-3 

i 


~ I 

1-3 ^ 


^ 5T M 

a 

a >-3 + 

w ^ ^ 

I M H 

^ I 
►O 

M in m 

*2SE 


ao? 


? 

? ! 

-* ** M *• 
M v | M 

9fe0a 

+ Z M + 

*^»M w 

►o a a »h 

^a < + 


2 o o w 

sg^af 

Z Z M ^^h3 

- ^ wl^j 1-3 V 

O 04 04 Z - - M 

oAAf 0^3 

4 o m 2 i 
os o o w wa; 
o m fn m b b o 

H Q Q a w 
K Q Q 2 a 


u y 


C-10 


ORIGINAL PAGE IS 

OF POOR QUALITY 


00 

00 

ON 


s 

s 

u 

£ 


M * 


-5 — « 


8 . SOI 


I ? 2 


H w. w 

y fa b 

O HI M 


ft 


l"9 

rH 

I 


W oo w 


S ^ 


a* 

* 


ft 


* 

M 

+ 




hS 

gs| 

VJ 


1S8 

?«» 

M 


it 


!2|m 


ro * + \ 
I *-»M W 

^*H > «->*ha > 

tigSg^sSs 

HI* M + H I * 

5 *-*~z M e 

I + b H CL I -b> 
*- Z "»*«'■» *H '- 
*n 'M . + hj + m 
*» M —»W -M — 
I w2 • H> M 


ggftg 
iSgV 
ft M • ? 

glSz 


zw*^«ca 5 i*^ >< 


8 


y « *•«. h 

b + b bT 

* b9 *.2 *> *H + 

M v M »b H + H 
wM — M «. w M — ’ 

HISslll 

+ * O + * 
M CN <n rlHWW 


M ^ J3 

as ^ ' *- 2 Z M 

1 g-ftf? 

| * 1 + " 15 '^ S ’ ft ' ft ' 
5 SSlSlASiJsMM 
¥¥¥“¥“?« 

*o b> os^j os h u o 

M M *-*M ^*5 
w wm ZHZhb: 

os os 


* w w 


ZgbhBh 


si] 

e! 


h > 1 

HH^ 

P 

in in 
X o o 


?dddd§8ii las 


~§ 

b> - 


'CL 

— £? 

b> -<n 0 


ft 


1-0 

iH 

+ 


~£ dT£?. 


h'W: 

srft £§ 

no *■ ' 

O * W M 

^hffln 

**88 

m m w w 


§ 

? ^rS 

■* w • • 

M M ft 

Hi 


^ w w 


m bi O 

a?*. - ~m s 

i s« *? _ 

S”I2^z^h8 m 8 

O w > b> wbj v? w qp 

ggitiptie ^ 8^8 ^ww 

M ? b^ W & W £ M M ~ w^M W w~- 

§ f?M S i£ifflH2H»bBbZZ2Z 

E '.+ fwfwo+w + p -w -S3SS 

EhhThI MOb>Mb^^MMMQQPQ 

Mores ' - - 

Bggfeo|:oaftjaa?H?-8-S 

| rr~§ ~§ if; ? I ?&sl s 3 3 

2 ? ? 8S SS m 5.i^ww ww 

u ^[oHHbb'- , aH(L'^fl»HafefeBh 

♦ MM wwwwZ ^ZCJS^SCJMHMM 

« s>fefehhQOQOyOQO 

♦ tit<HMMHQ wflwQwQw 


»-3 ha b3 ^»h bi •-♦ Si S S Si 


o y v o u 


y u o 


u o 


y y 


& 

S' 

04 


00 

09 

ON 



b> b> 


H 
s s 
« ♦ 

M M 

f? 

/% 


z 1 

w ui 

♦ ♦ 


i ? i> 



s 5 


g 

** 

igs 

z * ^ 


+ Z 8 


Z ^ X" 

nnStlB 

z s * z 

hRun U ^ 

z z m y g 

H H O 
H ^ ^ * • 

oo ^SS 

O O M • • 

SS 8 && 


* HnnH 


M 

10 

X 


ca 

ON 

VO 


C-ll 


OF POOR QUALITY 



& 


00 

03 

CT> 


ON 

O 


«-* ' * 
h) w « 
'Z * 
p~*hH w * 

ba wcu « 

p-* 'Cu c« « 

bgs? s 

hH 'P 9 * 

wCU p-*ha * 

cu Hh ' * 

pH 0* 'M * 

PISS s 

^ha CQ CQ * 
8 8 : 
: 

n a w w * 

tlis ; 

M M « * * 

W W • • * 

H WHH * 

ggi I : 

■ 1 ^ ^ * 

xp* #^» ha ^ * 

b> ^ ' - * 


~~ h> 


H M ~~ O *H h> 

'P 'p'pH pH M + ' 

W W W W « >b *H 


pH pH P 

X x <n 

p^hH P'hH 

pH Z pH Z W 

ssss S 

HHHH « 

6* fa fa 6* * 


p-* Z 

Szcaf 

S ' as as p' 

m * * ha 
ha — hj ^/*n^ v 

' -t- ' as - os ^ ^ 

H H H I H | - 

wZ M M pH 

os os os 2 os 2 

Ali~i~S88 

hahaOSpHOSpjQQ I 

' ' X £ I I C' 

M hH p'.hH p^R p-'P'ha 

-p'ppH Z h z b ^ ' 

os os •♦••'' w 

rrgSggssr 


K 2 


HHbbi 

H Ch Cm Cl* fa 



WHb + 

w + *» ha 

U T> M p-s p^ ' 

m ^ w i-a h h 

feHW - z - 

Em '■pZ hH - CQ 


pH w ha w 


s a . 


P-* ' * 


bg 


5 ~CU R wCL 

cu cs p*>» *>a(N 
hs '09 0 ,h 3 fcj €3. O 
*■ w p « ' Z Q * 

M SP f p-shH — « p-» 
w a. p-^ ^ wcu p-'h) 

a n b Z n h b ' 

« n O ' 'pH O 'hH 

p-s o* HHQ* HZ 


C3 S S 5 S 
b288bH'88 
SS££SE^ 

fflcpw wflfflyiw 

IlggllEg 


w P% W W 


Siek|?S§ 5 § 5 ¥fu¥?u 


o + R w as w p ^ 


p-sp>ha ha p-s p^ h^ ha 

ha h> ' tha ha ' t. 


hH P' Z Z Eh hH p-"— P-» 'PP'hH hH W wM M *-p w 

W b H 2 1 9 J_ 'P P' N S. *■» 'p- P»\ p^- '*»pH pH 'ppH pH 

p H+ 0 !HOSZfflHZHnbMbZZZZCdHHU 
E,HHftwhwQ + W + O 'W 'QQQQQQQy 

e= 5 <nw& si s Sbeb 2 .se s BasBSBBa 
ii Q il| e 2 e ^sfs^ara a s s a 
.AA gggg;7&;;3SaAa2?A?~g~l Js~B 


cti <T* o ‘p^hahsha p'pH b H SSSSmSSS 

SS S ? T S S S S hH hH hH W h3 ►o h^ hi M hH hH hH 
U Z 03 • * * • hH wwwIH w^-pww^wwwwww 
(0 55 h h b b '■p CU pH CU 'p CUpHCuCmCmEmCmCmOhOhOj 
• HwwwwZfnZOJHHUrgHHHHHHHH 


« HHwvwwZnZOJHHUCflh 

« ^(UhlUIUQOgOQOQO 

* hhHHHHQ'pQ'pP'-Q'P 


ps pn ha 


h> pH ha 
'+ ' 


OS OS OS 
H* + I 

p'. /sy 

h> ha OS 


S^hH' 

f Sf 
§«§ 
p-'Z 
X K S 


, , waaa 

Z 63 fa fin Cm b 

OS OS H W H H 


u u o 


PH U 


* 

* 

u u 


C-12 


ORIGINAL PAGE IS 
OF POOR QUALITY 


E-» W M 
X HI W 

o ^3 M 
ozzz 

5 88 * 


o^a 

j ^ • 

^gftsrS 

Hit*. 

assaa 

aseee 


ESftSi 


isass 


Iftftft SsSg 

i fa h b fa h O 


§?S 2 l§ 

«*s*H 

ddo din 

WXP WX 8 


> #— s-t 

s 


"S’ 3 

ii % 


■|J| |||| g 

**SjtelWi ; 

z 2 m g z z z « 

I iii'iii i 


ess 

Z I <o 

S p 4 m V £ H 

?rSH 8 ^ 

UZHD 
U W Q W -C 
WZhQh 
Z U '-m | E 
up * a h 2 

w r£rz%H ‘i 

oTSggSt 


: Tts ! 


S-SS 

hh z as 
? 6! STz 

Ila£ 

5 -~§ C3 

Q« ZS 
O * W 

?&&£ 

M X wfib 

rgSSi 

!§i.sa 

> M — -M 

> ou '-a z 
: x a x 

iH *"»«— l I 
I + HZ^ 

> M + M W 

> w(H Z Z 

> ft STaToT 

i xxx 


0 

X 

. 

X 

b 

B°° 


o 

os o o 

CL4 


• 00 00 ~ 

o 

z 

CSHHM 

Q 

M 

S? 

o \\§ 

• wnS 


C-13 


ORsGi^VL PAHc 

OF POOR QUAulTY 


00 

CN 


(fl 

0- 


00 

00 

© 


s 


© z - 

M - M 


*-% i 


i 

a 

0 




Hi 

2 * * 

WHH 

Si i 


eg:* 

sSsg 

H f?[j • 

¥tfl9 

M W 2i>< 

|i^ 

CQ H CL H 


CN 

cq 

CN 


R - 


£ 


m 

m 

rn 


r 6 


l .s % 

*~»© 64 * 

M | H _ 

&rij 8 u § 

‘ sras 

M U CL 


V * 


fi 

M 


I 


H M w £ 


h L 

§8 


S 


*H I 

!83 


64 fa g H 
M M Z M © 

fifi* 

63 63 


o, 5 g 


B ^ 

o S s: 

® ^ 4 “l © 

cq h i ao 

H H *• | ^ <*> 
w — M VO M '» Q 

64 64 5 ^ © VO £ 

63 63 Z on wS 

88EBS C, BS 
i 2 a iE 


64 64 < 


o 
o 
u © 


o 

© 


© 

m 

U oq 


g h[ 

w w 2 E 

*«s { 

64 B 




W 

5 


z 

w 

os 


63 


O 

S uuu 


£ 



i g x *a | s 

!@£ii - 

rH 03 m ^ i~ 


R 

Eh 

m 

8 


iS 

Siffe 

IIH 

&H W W ® fi 


oi a 

§ §S 

64 0* T. 
2n Oj 

IS® 


s 

3 ^ £ 

SB'S*. 

£«** 

CU X ^Eh 

m h 

S Eh 

e?5« 


* §15 «... 

M 64 04 04 — 'Z M &4 

64 ^ ^ O4 6 W j"H 

•* HCN CN *-> —- 

^^^^ 04 aCN 04 ^oEH^ 5 P 
OU O4 O4 CU ** ^ 04 v '•M ^ © 

- - - 'HO)aq*E<^^ ' 
fH oq rn a: ® ® m X E-» 

a £ aw ^T£i 


w «.«H 

1 m os 

[leg 


ip?** 

H Eh U j 

i^ll 

QWHHHM 


OO u 


SSSStK* 

:::::js*e? 

iiipi?tt§ge 

iSSSzzzzzSoa 

sssqqqqqwq3 


O 4 Q 4 Cu Cu 

® ® ® ® [ 

« « « « : 

as as os as 4 



c-ii+ 



o 

03 


<4-1 

ft 

M 

£ 

00 

00 

ON 


00 

8 


U 

£ 



C-15 




ORIGINAL PAGE' i 3 
OF POOR QUALITY 


g 

as 


2 

iH 

I 


ft 

fa 


oo 

00 

<T> 



V V 


_ eg 

o nn m 

O cn m 


in vo 
<n m 
m m 


& 

£ 


£ 

05 


eg 

in 

£ 



99 ^ s g 

as as < w fi 
- - “ M Z " 

- " M 

v v v " 

« * * g * I * 

£ £ £ S £ vo § 

t-H M M < M t-t 

as as os 2 os a os 
a a cu " cl q cl 


sis 

||| 

2 n q. 

^ « 2 
h 2 ^ 
b Q< d< 

-j cn ^-vr rg L> w 

K 9 c C? 9 c O* <0 

P P a* p -« 

3D "5 Ph - Ph ^ 

l CL CL CL CL '•"CL " " W "P> 

‘ " " " 'H CM a g* ^ — E-* " 

OHCMn 05 « "OC»hZme-i 
! 0 ?«¥?PpHp"-Q«H 


O fa 
O " 

in z 

AS 

*C M 

s ^ 

"«H eg rn 
o os os as ps 
>00000 


w © CO 00 CD 

I « « « « 

; os as os os os 



^ppppfafafafafaco^^ 

£EfiKK*Wi 3 »iB 5 

w „ *H CN fn ^ 

W M fa fa fa £h 

llllltlEBB 


8 


fa 

£ 

fa 


S 


leg 


2 as § S 

HHCL " 

"CQ D E-i 
W < CO M 
PdHX 

2 as o " 

M fa rf M 
-<o 2 05 

S 5 S? 

h m ' 5 * 

®.sse 

gg§“. 


*S 8 


gz 

v'Ba |85 



SSSSSSSm 

sVs^Hg’"" 

kbSe-zq 

8 8 8 Q EH 
fa fa W to D 


1 
eg E* 
to 05 

M - 

hS 
21 
rffi 
55 fa 

"CO 
to § 
*-> p 


!H 

I z 


fc 




||5 

w M w 
111 
f Site 

S SSfaB 

fa -"§ Ea § X 


3 ???ZaDCLfa 

sZ<co 




PS ""l 6) 

x in z 

258 
BS*. 

gs| 


fb e-i wo a: 

e 7 2 - E^. s 

M " D ^ 

■jmSS^ 

UEh^ ' i^o 

- M £h 5< H «* -* 
wP, 

B 3 ? 5«§2 


B -8 f ?B B 

Eh v. 


BB 2 - 

E^h Q 
W W '-'S 

ill? 

'S. 


BB£B *. 

._„bb's 6 b| 

B® rgoc.a-gilf'zS 


-P '\><2(NO 

^BSSSS 0 


*h cs fn in vo vo 



C-i6 



QR!G?t\ 7 !l P r )^T r:> 

OF POOR QtAuFY 


£* M D 

a HD *"* 


;?? -r?° i 

• 22 i A 221 o' 


I »* a + •*& 

ui c/j z «-a z 

X S 3 w 333 »h 

i J 1 1 l|£gg& g 

r¥ h¥ fwff&’f z 3 w 

a .- 2 2 ? 2222 . ???•???! 


► * W W Hr 


H I I 1 I U • • 

-a<-a»^ P P » 

^ ^ ITS T T 


) a 'Z vo a 'Z I I Z > 10 > OD + + + -f Zff' a «- O a s.mvoos 
jMMM® 5.® hhhhhohhcjhhojww 


r-t Q VO CTi ; 

<T> I W W • 

*H <! W Cfl | 


« I M M ^ 

<-A M W w a 

«-a a CO W M 
H M z Z — ' 

^ g S? ? 8 

fptn? 

Q ^Q5 M hj 

+ D * w a 
/NWAjgH 
J 1-3 « (-3 « w 

; * «. a» z 

-M «~*M i-H « 

' si S’’ » „. 

V I — '< Z M h 

'ggggz 


£ % 
8 ? B 

pLs 

g. g|g* 

3 §| 3 -g- 

BE 3 ~B*a 

a A |J) CM ^A S 

Zn Q , F 3 £h M /-» 

SB* iJSl^B 
ESS £§Sp 

a < 0 i a. cn —M 

*H 09 a^ 09 a> w 

illlsl^ 

a ^ 

^osos a a: w z m p« 
SBJppHS wS WM 
)p£Ea < EcuwHw 

* Ba 
rH 09 ^ \ V a 

l®UU«Mfeh<AAA 

* SS 2 £EbE-*B-'E-i 

e : os | < 3 < W A ^3 ?3 p; 

££S^.8 w 6 ^Eh 6^ 

S.SvS^^ v^H M M 


fi 

II I 

t 

*~s h 

*S B 

&4 Ckl w 

Ml 

Cm a Eh 

/as ^>4 a 

B^|B 

bb|§ 

ST &H 

S£B~ 

ttis 


Z Z Z W 2 

■ A A A Z 

|||g$ 

H H M 8 | 

%%%ll 


r- sep 

O A A M «*A 

O 

*■ HHhH 


3 H - 

ia.B'VUi 

iBHSpg 

! w j-wF, < 


15 -^8 PP^ h 
P^~ •■ «-P S' 


***« 


IBB^arg* 

!§g?|pj 

■ Is MS 


c!?s*5l? 

^ AO AAfe i 

a^aO £ &• ^ &4 
Em ^ F 3 F 3 am 

?g 9 BSS|| 
a® z § P E 


mm 


rH 09 m <+ IT) \o VO 


C- 1 T 



n p 

P © 


ORIGINAL ?«:•- i3 
OF POOR QUALITY 


S W - H 
DM X 

r ' “6* 


is asp 

- M Fh v« P> 

l«'S|8* 

ihoP 'm 
' £ x z Pin 

jSSSgf 


^ Ej § 

- m O rt e 

B *8SJ 

< -D 

as ■«• * 

OP6j£ ' 

bsb^ e 


H v[5ioO' 

.ag°^2 

‘g M .x£«< 
a \ -a 

<X fiCNOC 

as < 3 os o c 

W g wH «* < 

llllli 


gg|SE 

il!l 


; - ^cl os w cn * o • 

ip XQfHDZOOO 

!mos<Ssxd2omo 


• I 03 00 
. o • •* ««r 

O O O m m 


S Z S B Y B m t 

8 i8§888 g sSsfc£ 


* a d 

*-*0* M 

So 


(Ml H CL. CM 


< 

8 

• co , 

ro m i 

>s m m i 

M 5 M w- 
Pm h § ! 


d 88 

S . Q 

* ^ wo 

si? 1= 

H«Hn§S 

+ ffl wKJ - • 

m gq \\£> as 

M w . V£> ^ h 

hvo * a, 
P C9 I vo 5 

x as < -few. 


(«n irt so vo «— t 


W 

- 

« 

b « 

M 

3 

w 

1 

« 

« 

« 

* 

2 « 
5 « 

w * 

• 

§ 


ca < 

l l 

I P3 £*g g 

: cn z 
|o;mmm 


iSE i 

*bbLI 

ssslis 

O 5 u as h m 


Si 4c 2 

O O -o 
MM O 
W Wfflrl 


SSii 

-- B ll 
3 3 3 £ i 
33§!g; 


h o) n ^ 

os as as as 

£822 

CL CL CL CL 

CD CD 00 00 
« « « « 
as as as as 


»-k ^ 5 

SSp*. 

Pm 

{.n S8B^ 

d«;h a. x — E“ 

3|g 3-S- 

BE* ~S*S 

SS§> JS^P 
dgp K»p* 

h CL CL '-'Z m gH 
^ "CL CL W ^M 
M CM ^ eg -> — 

asBBSB^p^a 

O O M M CL fc v 

SSaCLCMCL^e^^M 
CL CL ^ ' CL - ' M - h 
' 'M (N CL ^ E-i ^Eh '• 


B 


b M Fh 
^ — M 
H CL M* 
M ^ ^ 2 

>tss£ 

— >H 

g*|B 

P Pp z 

M^ ?Q 


-* gs gc 'KhShp m^> 

p22£2STwm w 


s 

MCMm^r\^vv > 

illlsl^^ 

jllll^ss^ 

j Z Z Z Z Z U DCL STs 
QQQQQwSSwe 


WP '-'P 

2S£~ 

*“PB 

m^PP 


oyuuuuouoo 


C-18 



a. 

S' 

a. 



•riii 

j(N U U U 

!§§§§ 


£ 2 . 

*ht 

f WO 

sss 

§) &) O 
as a '.a 


o o 
«a in 

88 

88 

S 2 

3 S 

*&- 


_ *H 0> 

8 -lo 

8 1 I 

<*> m 

§1 


0*88 

*H X • • 

-. H yss-.gg 6 i . 
eSfgESSAtlfSg, 

<n u q a i 

MMgOMMMZMqS<Ot 


£ 

w 


b 
O Cd 


-H W b * 

V 8 8 1 ® 

d B U ncN 

StJir 

:n|li 

o < I ffi a i 
H I 6 < IM 

aiilli 



u u u 


o u u 


C-19 






I 

a. 


ao 

ao 

o> 


§ 

OJ 

3 



2 0 2 
' I 

Th iTd 

^ ‘HZZQSNOJZ <n T) J 

^ H “HUJ 'H hffl «.M 

sleissssat 



u uho u^o ynu y-H 


C-21 



x -pin 

v © m m 

DAS” 

« < « _ 

-iSS8 

« 


-?? oil 

M H f H H 

M < MM 

/> *-%. S D 2 

CC CD g a 

^£2 
nn o l x 
• • ti • ^ I 


H M I • ^ 'HHHHrlHH w M 

pM^5«M«^MM M ^ M M W w 

)hnfa£bMNbhHhb?§ 

hh^hPhhi<hm<mhEhS 


? ®is 

; sri 

• Eh • *—* O 

lf*81B 

• • + • M < 

I M M M •*■» Q£ 

( | W M W < 

rgAsrg '. 

i < M M fen M 


•hzs 

+ Cfl l *■'•{ 

* 5 o ! 

«*•« *-i Q H 
M O " • ‘ 

W M M CD C 

z S a? * 

| a es a : 


IS 8^ 
Dg*gt 


l§i llllai^ 

os « ^ -cuosw^m^ 

BRB.aassfa'?. 

§ IgaiUSE*^^. 

| z s.\p p | i A A j 


< <N <n •* T E 

- a. cu cu — ^ « 

• cu cu a «*• »^> a * h 

j w w *. £ m -■ 

■AAA” £ li l 

igasg^g M < 1 

4 a a a Qj » x~» r 

< S (S IK Q £ ^ MMM > 

• 2223 i • mjj * 

I *^0 | Oh Oh ^ " 

f "3 • <^^>2 SH > 

■ H CN n O ^&3 ^ Hj • * • 

iSggg?*j:££^^i 

ihhhbSbHHfahbr- 

IHMMM<HM<MMME 


I^.h'hI 


! M M < M 


<n O U O O 


C-22 


E« H ~Eh 

W H h H 


> E* E? h m < 

;e.VS| | 

jHHjg (0 

IlHe I! 
!??!& 1*1 
flH &Tm S' Eh § § 


lESB.' 


.\\ HH! 

: z z z 


« &H ™ 

Qt3^S 


III 

\.sSES 

lllll 


Pg tf ! 

« M ^2 
Z *’' > 

S*w 

to 7 p i 

7%^, 

I M 
S3 ^ H 

: 

1 ' ' Z 09 ' 
HCfl ^'-’U 
>Z wto w c 


> ft in in • 

) P ,• '+ fc 

( I + H •''‘rf 
4 ^ o** 25 ^ •* 

[nn g w 

^u)p« £ 

> 2 X 09 * - 

< >* * —a n i 

< win i • «* 

- | . /Nwh 

^Z^m 0O 

i P § s i ^ 

ixww>o 


*ys8 


-J 1 1 41 

Itij? i 

Z M M -MM 

islip 

vis 


09 wH V) E «*»*"» 

03 *^tO ^ § f* p 

*— * H* & - 4^ M ^ ^ 

•H(N£4 N • — M - 
09 X * 5 * wf> * W ^M 

*5 CN tO 1 O W w 

^ ^h h »£« S So 

MW- wf>9 55 ' I 5< OU 

' — M W v | W« I >5 

HWJwnH^wh ^1 

wp to £ wW to z M ~ 

I?8t8j8£ '-S 

55 ^ I !3 55 w+ to m - 
I £ ^ss + V3 <—6 wM 
09 '^pna, to w 
03 - - w m a, - x p to 
^ M M w H 4 (tO 

H w wtO W 4 w . X £ 

toppSvTojpw «* w 
z ~ TT •J'r'S'TE 2 ° o 

'll I E '~*5 '^*1 

H->hHZ ^hrtZni^ ^ 

I H Z ' '09 Z ' - I ' O) < 

I H - -MM - - M W >"3 W - W 

M M wwM M ww wM 5 I 

ICNwwtotO tO 09 tO w g ; 

ISggSiSsggSjgiSgiSg 

;totoxxtotoxx_xtozi 

is 2 8 


ifllii 


B 2 
ESP*. PS ts 
psrip 8 

*« tigs* s g . s 

§ Cfi W CU X wfH *T 

Qphi 'r 1 M «!*• ~ 

ci| Pe: I 

' ' j 4 H ' r* ^ n *H Eh 

M QC £ e ' M b &H Eh CU W 
iJ O ' £h w 'O m '^»Z 
m £ * M tO Eh ' 

fa cl a. wz M Eh to &H b S 

^'09 ft ^.sT W 

sgSsg^p^a 

)fc£cu£«^p3'^*' ^-*wX %. 

)§a09§^E^Sp 

« ' ' 0 < ' _' hh _' n .^ 5 '' 


jH'TtCLiCLCUQ* ''Q< '' hh ' 

Z Z < - - *» ' 'HC3 S^Ei ' 

OO^OHOtn^ag -OS m X w Eh 
mm oasasasasODHO—'B'-'M 
.tQWtDMggggS&Cb&atOrHw 


h M Eh 

Eh aT w 

M — Z 

tq Q 

2 vH ^ 


SbbSAsssb^^ *BS 

QZE«M H 030 ^\/- - - 

QhhU wOD ODOOOOHHMMfaEH*-.^^. 


JQQaS 4 4*4 
- - . asasccasos 


^•wwg g z b J<( 

lllil 8aMI 


5 8 M p P p 


to Eh r'fcH 

2^i;p 

W K wp 


izzzzzl a S’aT 


B P P P 2 ti 

w M M M \\E V) 


zzz w ; 

S sS S Pj 

5S8§| 

SSll! 


IH HWO ^invovo M 


<J> 

ON 

09 


u u u 


C-23 



K 

~ Eh i 

Eh M Eh E-< M 

Fa w m w w 


2 Eh ^ ( 

" 5 ^ i 

«-% Eh © ^ *» M * 

Eh m o Eh Eh y't* i 

b j^- *hT fa M 

§ P ^P§PP^5^ S | 

g S £ £SSeS£§ 8 

3 * s gr*s' 5 r|® @„- ! 

ojh^< Z os 

i *hES§PPS§^ - Q |j j 

! 0 & pP[t ^ P© ppPj p Eh M j 

isp^l^ooppsr^l g. 
r M 3 .Ssag§gtl^z|Sl 5 *1 

& EQ^OOnlD> 2 SZHK Eh 

m -faz 5 5 > a P 5 Si 

saBSSsalr^s las a 




I’i"’.’? D> 

?*v?b aa 

M M M M • • 

! p & > > § 


|P||Sa§||| fsfjBos!^ 


<0.£££DD>>U^ 

- SeSS 1 I T T X^ 

V<\ZZ ' 


os os - Kha^i-j^ 


ggSTTP^PSPi: 

!§!oo£ti££ 8 ?i 

9 9 Q°° D D> > -Jit 

88888 ^^^! 


*"* CQ CQ CQ CQ ( 

rr 8888! 


gfcgsS* 

D > * * * * f 

in in PPPP 
I’ f M m m w ■ 


^ gggg; x 


3 >4 >5 * X l 

i a > i J i J * 

4 in ub x S £ Q 

■N . . « « « « > 

HI I ^ 

3 ^ Hj ^ Hau 


I H • w w *Kh wwvKm b 

ia ^ar^sr r Sm^sssEmAm 

' ^ su • n in fa • w wm 03 o) w w w 
► > H M b) * ■ U H b D ^ . -~03 * 03 

:Slls^lls^g^^gs| 


rgsag"”: 

- O 03 0) Z Hf* nr : 

IfaSTfagSSj 

I M M w z : 

888 


! P z 3 2 


i?uhAI 


!S rtH ST| 


3 2 2 

m CLi CL CL 

'O h n n 
O OS 05 OS 

>rggg 

]Haaa 
: m 

] '—•00 CD 00 
t * « « 
os os os os 


H r* 

5 * 

Be s 

PPpS 

bHh ' 

g PSrH CL X ~&H 
QO iJ m -'•-t M 

si | 2e a ~s 

2SS 5 s ^e 

wgV 

fa CL d ~Z M Eh 
" -CL CL 03 
HCN <N ->—* 

^ i i ^ i 

SfafaCQfa^P^P 
cl " "Cl ^ ^ m " b 

-Hd CL **• Eh -E - 
hTBSCC ^hSh^ 


KOOHO wS 

SWW^BS 

^ .H 03 cn «<r ^ v *. 

I 11 111 

EEEEESRhVe: 


P 5 ti ST 

ft s 

51 l 
rr * 

PST w 

M - r Z 

2^*:- 
»te tie 

?b|e 


g03 Eh ^Eh 

tlR^ 

“ B ‘S 

Z Z g ^ z 


Q * OS H 

2 Z H X 

M M fa " 

^ W 5 Eh 

832B 

8*i«H 

z as u - 

H CtJ < M 

S 5 m P 

Ills' 


SK^s'SE 


SSSEESSm' 

dV«**mS 88 

Hill''' 


QQQhSjZ* 
OQQ OS W * 


HC3 HHHH HdH 


C -24 



* 


id 

o* 


oo 

QO 

a* 




i||l 

dll 


.afiiS 



w 


s sss 


k 00 * 


£> W QO 
w m 2 

ass 


? 

& 



\ 


« w 


8 


Id 

VJ 


R 

B 


a* 

03 


1 


ff BB s 

MpE^S 

f-OMH . 

«f?s 

£h £h v ^ 

b Eh 
« CL *- -b 

-SSSe: 
fi -:■!£& 

-rsp s *a 

«■* *> **■■* «*. 

BSK^fe 

f^HtH e7 

M a U wM 

ll^Sl £o 

iilgariill 

QO *• - - ^2ZM< 
cq <n «h 


OOOOOOOO 


o 

c« 


o o o o o 


n 

£ 


S> 

S' 

CL 


CD 

00 


00 

IT) 

s 

<# 

«N 

1 




R 

E* 


O 

O 


CO , 1 . . -, i 

^ © « 2 a a j 

g^j^^oCLCLCucu 

O O M cn rn -* 
MM O 05 05 05 05 
WWODHOppO 
22* J SlEKfi 
PP^mcucuclcu 

5 M S «ffl 03 00 00 

Q Q 05 * * * * 

- - 05 05 05 05 05 

ssslii 

M M M CL U U 


V o u 


; M Ph 
v CL flu W wm 
oj cn 

sssRwa 

W A & A •"> b <"** <*> - - - 

SS^g^E-'^p S>|B 

0505 ^ 05 M Z M Eh 

P P M p ~E wM 

SSSJSSaa 

lllll^ge 

lilliigtg 

z S 5 a ST 

Qw53w 




r 

B* 


Eh' 


R 

Eh 


R 

Eh 


r- rr? 


Eh Eh ~Eh 

M M Ct< M 


n 

, K 

Eh — 


Z S I 

Jpg! 

£> h^O |h H*£-> £* 


I 


jrSRRtJi 

w ; t E* X 

M O < &H Eh M Q 

gSg^gg 

BHBPnVPSRRJiB ! h h 

& "e^Q w«o® Igf^ZH^v 
m - CL Z 'r 'ZZOHpZ 

S '— H^ZM-sEn— '<52\ X\Q 
Z tH M -><E^\\a Sb? 
M u b O H"" < \ - 03 ff>M < < b 

-U v\x2(NOQpn 05 05 ' 

~ §^S 


HCNn nr U*» Vtf wo 



C-25 



£H * 

s * 

ffl 

8? B 


^3 CL CL CL CL 
rf ' ' ' ' ' 
— O HCNn ^ 
O 05 OS OS 05 


3 OSH Q-X^E^ 

Q O H cn 'H M 

sij ~b°J 

fa a |al? 

h E *• HWE " 
b Cl CL ^ZhE^ 

' 'Cl. CL ffl wm 
H« '•*«• CN -> « 

' OS OS 0£ o '-‘U ffl 
PpCLp ' Eh '*C 
GL GL CL CL *^*0 **• 

5 5 -5 ^ -e* ^ 

D-CLCSCL^EnbEH 
' "CL - -m 

-fH CN CL E-< ^ 8-* v 

OSOS '95HSHE 

2 p h p wm 

£ n £ it/iH ^ 

WAS»i8 § 

H 05 ro «<r ^ ' - 

Himgftg 

I5I5IEbbb 

zzzzzl 5’S’oT 
QQQQQwSSw 


w 

R* S 

51 e 

BB 6? 


£ 'H 

>tS0B 

b^|b' 

b||| 

w e* ^8 

gog,; 

*°gS 

ills 


; m '•e-' 

ix h < 


lx C9 S 


i* 26> 

I O w z 


““eats 8. 

s > a ,'S < & i 

o, ftj a ph z Q m ; 

1 1 ! Esses! 

iiiiiPJ 

ffl ffl ffl 8 fe 

||||| 


< EH ~ 

CJ M 

ffl wE-» 05 "■» 
• ^HbOE 

l Eh -M Ph X 

!2E*H '2 
-M b M Eh 3 


!&B«aS°^ 

l M — vZ M 


E* vO ^ Eh Eh 
h '■'O Fh Eh b M 

x x^X 
M O < Eh Eh M Q 

r? « m m x s 

P~®dST|° 

xD ^ W tO < ^ 
< < -Eh 

• Eh 2 ' x -v y-x H) 

t o xo <^^«E ' 
.1^0 EH Eh fS Eh 

lEnf^hh 'H 
M vw V vE - 
-o ffl Eh Eh H Ql F 
I jog M M W 2 & 

foS s P > 1 F £ 




ZZEhm- 

h U n O ^ 
- u -<x 

VEh « 2 s 
Q M ffl 2 


^ S \ 
X 2 *«• \\0 


:^UZZZ 


i Eh M 3 Eh ' ' 

|||| ijplg 

iii ill! 

Sis.allll 

si|BiiS|S 553 

s ^ ^ § s U fi 4* A A 

Q ^OSQQQQClclcl 
0 £QQ0oooaQ.a 


*h U"> V£> »X> H HCfl H 


»H »H HI *H H 


£h W * Eh 

&5 S ! 65 

?t!¥m 

o z * h < * 
o F * o i i 


: Bz V 

1 -3 £T 

f f§| +' 

s — * - ^ 

t o * Eh M 

r ^:s £ 

2 M * Eh * 

^53 w h ^ 

J M ffl * H H H 

'Eh + * I I H- 
) « H H < 

M * H E | 
wM w* K 
■ * Whg 

| rf; W* CQ H vZ 
I I « O + H M 
H M wL 
! H H * Z H z 


'OW 'Eh 
tr* -z Eh 
Eh m - 
b Eh X w£h 

:«BB°*a 


SlT* S^ffluQ^rtuS 8i£8 


H fr3 - 

i if 

IglJ 

M O '—O 
y*’ M O 
> W®H 

H H ffl w 

J ?“a 

g s &e 

Hgu| 

z B j 3 § 
g sftfta 
gig- ,aa ' 

a 63 w 


^6E<P« 

— OJ l-l b &• 
'•^^Z M Eh 

<h£ q ^ 

bsbJb 

Eh <S W Eh 

M Bl IS 

-MU ” D Z 
£8-^ 2 - 

gzjBBBV' 

H O Eh ' ' 'Hi 

ffl £ w Eh Eh Eh Z 
Z S ^H M M | 

eoE 5 o*cl 2 
h wsawH 

Q ' ' ' 'Z 
h esn rl 


I (fl *— * 

I Jg lo »-3 ^ 

H H H M 

!: laaa 

• « I I I I 


► ffl ~Z i 

'^§^zi 


I « < ffl U Q 


g N.^ jjO P g 

3 jj T TJ oJwhk 

:g<ul84£8i 


ill 

§g HH 

&c vi 


C-2 6 


U-) W 


09 

00 


00 

m 

s 



O fa 

S 2 
^ 3 


N ^ «— I «N f 
I O Q 5 Of { 

o'fpps 


- =£555 

*h - T r> cu a* a 

§§^8 its a 

gSS?X2£S 

tf 0,01 

p H H U —00 00 09 

aLLsf ! 

M g g U i 

o plis 

u w 


uuu 



S' 

a 


09 

09 


09 

m 

8 



Hfgn 
0 U fa fa 

SJttt 

US US 05 <n <* 
,«**» *-■» cu cu cu 
PO C 9 CL 0* fa 
fa CU **-» *-» 
fa fa ' 



• • a a o. ^ a • 

*^444" S£ 

^ . . H M H o W < 
n rfl h ^ h 
CU CL ft «4 A < <*—» 

“^“rsiss S3d 

CO CO CO | 0,0, 

• • -HcNno^b • • 

1“ ‘ SSSSftid^v 

£fu£««>S«iE E'C M M 

bhhbtUIUIUHhlUbH 


fa »H fa 

H e-l M 



C-2T 



I 

cu 


03 

00 

a> 


oo 

in 

s 


I 



I 

cu 


00 

00 


u 

i 


2 <2 & 
3 p ' 

h 

Pu CU CU 
vCU 


* 

Eh" 


* 

6 "' 

M 

o 

8 e , 


S 

f<‘ 


SI (ft . 

9 « rH CU X wfri wj-s 

3 2 a, 3 ^ > w 

sil ^e®gg 


2 

cu 

00 c 

* r 

OS t 


M W E-i J* 

- STw 

sBsss^p^a 

) g g cu g £*»d ^ 

>a.cucNcuhHbp 

* *> ■>‘0. •*> «. M *.Fj 

- Cfl CU E-< „*• 

Ci ' B! H J H 
OpHO 

E&cuEcuto^^-’ 

S^^ 3*85 

. . wm 

itlltt-SSB 

i???S 3 j!g££ 




C-28 


£ 


03 

00 

a\ 



iiiii 



M Zb 
h Z M Q 

*• H ® h 3 m <^»<^*»H 03 *H *~»*~%<*%*~»*~»t~**~s*~%**s^<.*+*^+ 

O' P J MHHCQNNCNCLm^»tr»voc^a3<n<®yQ 

is _., 

g? s r?n gggggggegggggggg»Ai£?s 5^5555555513 

^.Ljf , ‘' , '' , “ v ' s '''' ,, ‘ vv *‘' , '‘'“' > > 5 l x-'p^r-'.f-'-r^f'^r^p^r'^r^t^-C'^r^ 

v w Qa>a\asa\a\asa>a\asasasasa>asa\<T\a>asa\ w z z — - w q wwwwww.wwwwww 

g 2 1 q g a g 9 S S a 5 S a S S S S S 5 S || A A g S|ggg2ggjg|!i5g3g2g!g 

8siiii!ii!li!!l!lll!il§il8§iiiiiiiiisiii 


MCNftn' 

gggg! 
uuuuuul 
^ >n <-> ^ ^ b ♦ 


gg 


u u u u 


u u 


I 

£ 


® 

00 

as 


8 

*• 

CN 




6 ) 6 ) 6 ) 65 63 6 ) 636363656)65 6) 65 63 



o o 


u o u u u o o 


C-29 



Tdl' 

' H H U) Z ' 

r s' § a 

rsgssf 

isasis 


3 ^gigg~ &a ; 

4 Htnirunn J mm h 

* ^ m ». ^ cl &*&•*■"'"■ 

:oa<M -Mm£*-* -ZPO< 

jygAVAiMSSBJS 

JS'-'^^AfLiSS' 

iZ^fO <n «■ h h t/i Z / 
pea ^ ^ •<-'» Z Z - 5r 
wca m ^HHnSS l 
J Q — ^ Q - 

I 5 ^Sd|S^SzS 1 I 

!Sg££g££S£gBi 


s » ^ I - «. *— » >*. 

> --'»© O a 1-3 *-> O *-* 

lOinin ^ 3 3 05 in o 

unnn J mm vnin 
*. cl b f ^m 
i *»M <-5 £ ^ ^ZM^m v 
Mil v^. M M M «C 8 M 
4 I MMMMMPJ km i 

*• M -^IJJQEhX -m 
■n ^ ^ m eh a m w x-* ^ 

4 ^->«H M ' M M (0 2 oa -* 

" M ' -‘'-'22- 5-^* "C« 

H ** M M M 5 5 I — 'M M - 
■* m w — -.wwcnQriwM 

!? 8 &£ 8 »& 3 ga 5 ’ 

U&iiSSami 


« a e d 

3E V £ 


a 

M 

§ ff 


8 P 


h m o p : 

siiv 

iiigi 

m cj z i : 


>4 f 

‘S 3 

IIImm 

ijifi 


^ i 

'QU - ~m a 

>03 m vo - 3 

M w . M M 

m w a 

h ^ in ce m EH ■» 

•t J ni Q X m m 

> £: A m m o A 

4 • M Z I in x -~»M M 

l £ 2 D Hffl pin « z z 

O ^Qm^H§eT^ww 
fc ZiawMrfwH 
HZHCDQg£Q(/]Z 

ST“ SoSposSS 

mODSQD5O£05UO 


m a * 

H fj 'H *- <- 

® g ni 


- & ^O rf 

i <--* ^oa m c 

; m oa I m ; 

'i J f -- 


I ~ — Q vo c 

-a z in - 

S i w Z m m ( 

MS*" 
a a Q t 
i u O os 0 c 


J M h3 M 

^ ^ a % 8 

X 1 H JL 

m a ' *r 3 < 

"3 M " MM 

M £♦ ^ r M fa fc 

H ^ '-mm2 M *-*■*• 63 

inn j ns X m cn 5 

-AAB^ ~i£S 

«-''M M Z I ^ AMMtfl 

!\§§t w ®g\g§-, 

:ra — w Q r-. m z m •— 

iQWZMMQMQ^Za 

ilssgoiiiagb 

IKU0KQK0CU04 


c-?o 


§ 


25 *. 

5 sails s r 

r* M £ X m 



S 3 


v v «*% DO U Q3 Qe 3 Cd 
^ ^ ^ - H -MHOUijaU 

f<i— I iHJEJCM I © | + I fH 

§ M i » C C ^ w i ha I ■ I . < 

I P o> h w h *^ ■»■'»'»»•*. 2 J | | / 

'-’ON fllbh^Z^ZaiZHHHHMHCT «-*l- 

$SSmm^E^M^M^DD>>£:S 

s s a g s srs s ?a a %£ a 1 8 a =rs £ s 


ON 

ON ON 

uuunuuu ooouooG 


8 

d 

M 

B 


on ao * « 

oyy?Suyyuuuo 


— *. « 


W W 
»"0 

- - WWb ' 

^ t-j «»fH 
H WH Wrl ^ *- *H W 
tflHWHWHHWH 
*3 "b W W H - 

siasspll 



D v l>'§c:ggSa 


£ H2H2 2222 


OOQSMQSmDSQSQSQSQS 
incuftsaascuacLacu 
ON _ flu _ CL 




uguguuyyo 

O ^ 

u^a^alii^a^s! 

Zmmmmmmmmm 

afahhbfahhhbb 

PmmmmmmmmmE 


8 

8 

M 

e* 

M 

z 


I 3g3g33333 g 


OS 

a* 


& 


Isis 


GO 

O'UUUUUUUU 


g> 

S' 

CL 


CO 

GO 

ON 


® 

U'» 


CS 


O 

ID 

ON 


m x-n - v a. 

?T*h *f T 1 ® m 


O W J — 'O ^o o u ^ 

in «-* 3 < m in o in tn -> 33 ^uioi 

m t-3 mm xninnn^ mm -mint 

- * 

M * <? *? £ J 5” |H M rt pjrt 5 Hlw 

~A A Ap£~ M .~~A A A8E*~"~ 

M CN § § “| 

w 


o 

in *-* 

*n J 


•g^T. 

wCN M 


^nn xhh W2 


M M ON 


SU M W ON 
“ woo 2 Cd 


Is! 


^2 I 

d r ■ 


sss||§ips§^s|i|| 

o < os OS v/^^ww(J 04 !BS 


<n§§ I w-nH^HH?? 

HM^S § 


n wm - 

Q%Sl 

pi S 


a. m 63 
u ; 

1111 


HHHHH 



ao u o <j u u u 


u o u 


C-31 



CO 

03 

<T> 


S 


10 Vi 

- 10 Vi 1-5 ** 

<H ^ ^ '• *— I 
05 ^ ^r-» W 

HHH 71 H 

- 10 10 M 


gggSS 
3 5522 

S HHHH 

z z z z 


« 



0 **>»*V*i fi * 


P:\CE iS 


OF POOR QUALITY 


E 


* 


u u 


-o 

03 in 

in cn 

m v 

" o 
oq in 
in m 
cn - 

EH in 


I 


§ 




a.- - . 

8 S T 4 8 o | 


o o 
03 in 
m m 


—.CN 
tO 

«n wb 

tO CO - 
-03 ^ v(N 

i to 03 - -m to 

) H to IN OJ to H 
co to M *> 




»BS83 

gggggg 

mSSSSS 

Z 


fcggggg 

Z M M M W W 
►H os OS OS OS OS 
OS CL CL CL 0- CL 
CL 




H 

a 



3 

s 


uuu 


S' 

<0 

CL 


03 

CD 

C\ 


s 

s 


o 

o 


se 

§«S 

Tgii 

813 



•-> »n ~ 
hhS 
Z Z m 

ss 2 . 


in in 

* « 
« « 

M M 

z z 


CL • CLLUO HhO»-5 

£|£T££ M ^S1 

CLHCLMHZlOlO I M 

m os I I m m a 

d 5.j ^ 2 2 ^ i 1 

2 &T< waSQQPP 

UMUMMHjQQi^ 


Hh /N« «k M ^ 

z z *Z ^ k '"L a 

hhhhcohhhhU 
1 1 1 1 wwqs 

2 C? Cv *** CL ^ «—» gj 

h m m ^ *3 nmij^ 
=-•03 00 ** ^ »4 CD oo - v 

388 > du 88 >du 




UfabhUOU U 


> « o 

* o* 

u O U 01 


a> * * 

Suuuuucj 


c-32 



CALL THRUST 

IF(INCALU) CALL PRINfr<l / l / NI / NJ / PDUHMY,HEDS, IS2, JS2) 
IF(DATAfc.BQ.BTEMP) GO TO 56 
OPEN ( UNIT-53 ,FILE-FILOUT,FORM«* UNFORMATTED* , 

ACCESS-* SEQUENTIAL* 


OJ 

u 


£ 


^ M 

£ 

X XX 


M D M 

w M | 

o *** *» 

-.1-3 


£ £ 


h rr, 


~r£3 * 


ST A 

** 


c a T* T 1 fe 
5 t'tr 

^ ^ n » rH s. 

M t3 JT| | Hh 

-wai 5«as 

"*■ *— I I-H *-» O w 
^ I I H H J 5 rl rl 

S bb I 2 EqiJ 
* “O «* 3 —i-3 >*3 

“ (— I %. .. 

1 2 25 mm x^^mm^ 
rn m s. ^ Z w m 2 2 r 

in *h m v «*. 22 *“ •‘in 


-mmbi-3 hh | | i £ 

-» M “ I I M M M < 

- k M MM^v2 Ed 

O C3 v s. w \ *— s fij Q£ 

S '-'in W >"3 >-3 *-«•>* CO ^ hj »-} — * 5 M 
B«HH vgw^wMMMM 

.lilllllllHpli; 

I MHHHH HMMM 


I 

i 

i 



•W53S 3 

'»>D Ol — — --< w 

3wx||g'og| 

_ X 2 2 p § 2 § 

o £ 2 ca 

l «H 09 CO ^ rH 

o o m in 

«f M MOO O 

in oin ro ««j* tj* -*r 


V u u 


C-33 



